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0.1  Objectives 

Specific  objectives  of  the  research  at  United  Technologies  Research  Center  (UTRC)  are  as  follows: 

(1)  Develop  methods  and  tools  for  uncertainty  management  in  non-linear  non-equilibrium  models  of 
wave  phenomena  in  the  jet  engines. 

(2)  Develop  methods  and  tools  for  parameter  identification  and  model  validation  techniques  for  non¬ 
linear  non-equilibrium  models  of  wave  phenomena  in  the  jet  engines. 

(3)  Develop  methods  and  tools  for  Design  of  Beneficial  Wave  Dynamics  for  jet  engine  life  and  oper¬ 
ability  enhancement. 

0.2  Summary  of  Accomplishments 

The  accomplishments  of  this  research  program  emphasized  the  more  promising  directions  of  points 
(2)  and  (3)  listed  above,  where  the  uncertainty  management  called  out  in  point  (1)  was  implicitly  ac¬ 
commodated  in  the  analysis.  A  particular  emphasis  of  this  research  was  the  identification  of  nonlinear 
models  based  on  high-speed  video  of  flame  dynamics.  This  is  a  key  component  of  the  thermo-acoustic 
feedback  model  studied  extensively  in  previous  AFOSR-sponsored  research  programs  at  UTRC.  These 
results  are  presented  in  section  1.1. 

These  concepts  were  extended  to  nonlinear  model  reduction  based  on  tangent  space  approxima¬ 
tions.  Here  local  gramians  are  empirically  computed  based  on  perturbation  trajectories.  Key  com¬ 
ponents  of  this  research  were  the  development  and  application  of  algorithms  for  approximating  the 
nonlinear  manifold  for  which  a  data  set  belongs,  and  identification  of  the  nonlinear  dynamics  on  the 
particular  manifold.  This  lead  to  the  concept  of  a  hybrid  (switching)  locally  affine  dynamic  texture 
model,  based  on  the  local  tangent  space  approximation  of  the  manifold.  These  results  are  presented 
in  section  1.2. 
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Stability  Analysis  of  Systems  with  Symmetry-Breaking  is  presented  in  section  1.3.  Here,  sufficient 
conditions  are  established  for  a  nonlinear  PDE  model  of  thermo-acoustics  with  wave-speed  mistuning. 
In  investigating  a  new  approach  to  stability  analysis  of  the  nonlinear  thermo-acoustic  model  we  de¬ 
veloped  the  concept  of  continued  fraction  convergence  as  a  condition  for  stability.  Initial  results  are 
presented  for  a  string  of  oscillators  and  the  concepts  are  extended  to  the  application  of  subsystems 
connected  in  a  ring  structure. 

The  results  of  the  current  research  are  summarized  in  7  journal  papers  (1  published,  1  accepted,  5 
in  preparation),  9  conference  papers,  and  another  2  conference  papers  currently  in  preparation. 
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Chapter  1 

Summary  of  research  results 


1.1  Empirical  Modeling 

In  this  section  we  discuss  results  on  empirical  modeling.  The  first  set  of  results  relate  the  standard  lin¬ 
ear  dynamic  texture  modeling  framework  to  the  recently  introduced  concept  of  dynamic  mode  decom¬ 
position.  These  results  are  built  on  principal  component  analysis  of  the  data  set.  The  next  two  results 
focus  on  the  identification  of  nonlinear  systems  and  the  analysis  of  metastability  with  the  analysis  of 
the  Markov  operator  governing  the  transport  of  distributions  on  the  phase  space. 

1.1.1  Dynamic  Mode  Decomposition 

Through  the  singular  value  decomposition,  we  show  that  the  eigenmodes  computed  from  dynamic 
mode  decomposition  are  the  same  as  those  computed  from  a  linear  dynamic  texture  model  computed 
from  principal  component  modes.  Comparison  results  are  presented  based  on 

Dynamic  Mode  Decomposition  (DMD)  is  a  method  that  allows  the  extraction  of  dynamically 
relevant  flow  features  from  data  [66,  69,  68].  The  relation  of  DMD  to  eigenmodes  of  the  Koopman 
operator  appear  in  [63]. 

Consider  an  ordered  sequence  of  (vectorized)  data  snapshots  x*  £  RXI  where  k  =  l,...,iV  For 
integers  jf,  k  where  j  <  fc,  let  =  [xj, . . . ,  x^]  so 

X:=X?  =  lXl,...,xN]~  (1.1) 

The  row  vectors  z/  €  RN  can  be  considered  as  records  of  scalar  time-series  data  at  a  single  spatial 
location  indexed  by  l.  In  particular 


'  Z!(I  :N-1)  ' 

■  z\  (2  :  N)  ' 

Xf-J  = 

.  zm(1-.N-1)  _ 

— 

,  ““ 

.  zm(2  :  N)  _ 

Suppose  that  there  is  a  linear  model  A  such  that 

X^  =  AX^-\  (1.2) 


4 


Suppose  also  that  we  can  write 
where  S  is  the  companion  matrix 


■-[ 


In-i 


:=C  +  J3 


where  s  is  the  least  squares  approximation  of  xn  =  X^  ls,  so 

i-i 


and 


Write  the  eigen-decomposition  of  S  as 


[xf-^xf-1]  Xf-lTxw, 

,  cf  =  (l,0 . 0], 


0  1  ■ 

1  „ 

o 

II 

/jV-2  0 

|,  D:- 

0  s_ei 

sw  =  wn, 


(1.3) 

(1.4) 

(1.5) 

(1.6) 

(1.7) 


where  the  columns  of  W  are  the  eigenvectors  and  ft  is  a  diagonal  matrix  consisting  of  the  eigenvalues. 
Then 

>iXf-1W'  =  Xfr-1SW'  =  Xf-1W'ft  (1.8) 

which  implies  that  X  ;v  ~ 1  W  approximate  some  of  the  eigenvectors  and  ft  approximates  some  of  the 
eigenvalues  of  A.  We  refer  to  these  eigenvectors  as  DMD  modes,  denoted  by  3>,  with 


$  =  xf^w. 


(1.9) 


Note  that  the  DMD  modes  are  composed  of  the  normalized  columns  of  W  scaled  by  the  data  X^-1. 
Therefore,  each  column  of  $  will  be  scaled  according  to  the  data.  This  is  analogous  to  taking  the 
Fourier  transform  of  a  time-series  record,  where  associated  with  each  frequency  is  a  magnitude  re¬ 
sponse.  This  is  discussed  further  in  the  next  section.  Therefore,  associated  with  each  DMD  eigenvalue 
will  be  a  magnitude  given  by  the  norm  of  associated  column  of  <E>. 


Principal  Component  Analysis 

Principal  component  analysis  (PCA)  is  also  known  as  proper  orthogonal  decomposition  (POD)  [35], 
Kaxhunen-Loeve  Decomposition  [72],  and  singular  value  decomposition  (SVD)  [36].  The  SVD  of  a 
matrix  X^”1  is 

xf-1  =  UZVT,  (1.10) 

where  U  and  V  are  unitary,  UTU  =  I  and  VTV  =  I.  A  popular  method  of  computing  POD  is  the 
method  of  snapshots.  First  form  the  correlation  matrix 

K  =  xf-^xf*1  =  V'EUtUZVt  =  VAVt,  (1.11) 
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where  A  :=  E2.  Through  the  svd  relation,  the  POD  modes  U  are  computed  from  V  as  follows: 

U  =  X?~lVA~2, 

(1.12) 

which  is  identical  to  (1.10),  however  K  is  of  smaller  dimension  than  X,N-1  and  therefore  computing 

V,  A  through  (1.11)  is  sometimes  easier  than  computing  (1.10)  directly. 

Combining  (1.5)  with  (1.10)  results  in 

s=[xf-l7’xf-1]_1Xf-lTx^ 

=  VA-1VtVZUtxn 
=  VX  ~lUTxN. 

(1.13) 

From  (1.12), 

xN  =  Xf^s  =  X^VTT^xat  =  UT.VtVT,~1Utxn  =  UUTxN, 

(1.14) 

which  implies  that  the  least  squares  approximation  of  x\  is  precisely  the  projection  of  x/v  on 
modes  U. 

to  the 

Conditioning 

The  algorithm  for  computing  the  eigenvalues  of  S  given  by  (1*4)  is  ill-conditioned.  Here  we  discuss  the 
use  of  SVD  as  described  in  [67]  to  result  in  a  more  robust  algorithm.  We  essentially  apply  equations 
(1.7,  1.8,  1.10)  and  apply  a  similarity  transform  on  S  so  that 

s  ■.=  si^sysr1. 

(1.15) 

We  first  combine  (1.8,  1.10)  to  get 

auevtw  =  ut.vtsw  =  ui,vTwn. 

(1.16) 

Pre-multiplying  (1.16)  by  UT  and  applying  (1.15)  results  in 

utauy:vtw  =  zvTsw  =  e  vTwn 

(1.17) 

=  7SVTsvErlTVTw  =  s  vTwn 

=  SY  =  YQ, 

(1.18) 

where  Y  :=  E VTW  is  the  matrix  whose  columns  are  the  eigenvectors  of  S.  Referring  back  to  the 
discussion  in  [67],  note  that  according  to  (1.17,  1.18),  UT AU  =  S,  and  the  DMD  modes  are 

$  =  UY  =  UYVtW. 
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(1.19) 

It  remains  to  be  shown  how  to  compute  5  directly  from  the  data.  Define  ©f  =  U1  X,v  1  = 
EVT  =  [0i, ... , 0/^_i]  and  similarly  define  0*  =  [0j, . . . ,  0*]  .  First,  swap  the  left  and  right  sides  of 
(1.3)  and  pre-multiply  by  U7  and  post-multiply  by  FE-1  to  result  in,  after  applying  (1.10), 

5  =  EVtSV£-1  =  l^X^FE"1  (1.20) 

=  UT  [c/EF2A'-lT,xN]  FE-1 

=  fe-1 

=  ©^©f-^A"1.  (1.21) 

The  above  equations  involve  computing  the  SVD  of  the  data  X  and  projecting  the  data  onto  the  PCA 
modes  U. 

Truncation 

The  scaling  by  A-1  in  (1.21)  motivates  truncating  the  model  because  A  often  features  singular  values 
with  very  low  magnitude  and  the  high-indexed  entries  of  ©f -1.  In  other  words,  the  projection  of  the 
data  on  to  the  higher  order  PCA  modes  results  in  time-series  data  with  very  small  absolute  value, 
which  is  the  case  when  the  number  of  snapshots  N  is  smaller  than  the  spatial  dimension  of  the  data 
M,  i.e.  N  <  M.  Suppose  that  instead  of  computing  the  DMD  on  the  snapshot  data,  we  first  project 
on  to  the  low-dimensional  space  spanned  by  a  truncated  set  of  PCA  modes  Ur.  Equations  (1.2, 1.3) 
become 

©2V  =  U7X$  =  Uj  AXf-1  =  UjAUrQ f-1  :=  Ar©^_1.  (1.22) 

©^  =  U7  X%  =  Uj  X^~xS  =  ©^-1S,  (1.23) 

which  results  in 

Ar©*"1  =  ©f'-xS.  (1.24) 

Rearranging  (1-24)  yields 

Ar  =  ©f^SFE-1  =  EFr.SFE-1  =  S.  (1.25) 

In  particular,  by  (1.13) 

S  =  [©f-1,©f-1s]  FE-1 

=  [©^-1,evtvt:"1c/t*^]  fe_1 

=  ©2f-1,0Af]  FE-1 

•=  ©^  ©^“lTA_1,  (1.26) 

which  resembles  (1.21).  Furthermore,  by  (1.14), 

9n  =  UjxN  =  U7UEV7s  =  ©f"1*  =  0?-WE -'Utxn  =  ©f^FE"^*.  (1.27) 
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Dynamic  Texture  Model  with  PC  A 

Dynamic  texture  models  are  essentially  linear  dynamic  models  describing  the  dynamics  of  PCA  coef¬ 
ficients  under  the  assumption  that  the  driving  noise  is  second  order  Gaussian  (see  the  work  of  Soatto 
et  ah  [23]).  Suppose  the  data  snapshots  {zk}k**i,...,N  are  the  output  of  a  linear  dynamical  system 

=  AOk .  +  Bvk ,  (1.28) 

xk  =  <&9k  +  wk,  (1.29) 

where  =  [\^i, . . . ,  *&#  ]  €  RnxK  chosen  to  minimize  the  objective 

(1.30) 


Note  that  in  [23]  the  minimization  is  based  on  all  A  data  snapshots  and  here  we  leave  out  the  last 
snapshot.  Under  the  assumption  of  stationarity,  the  optimal  solutions  will  be  similar  in  both  cases. 
The  minimizing  solution  of  (1.30)  is  obtained  through  the  SVD  of  Xf^1  =  U'LV1  by  taking 


^  =  [/,  ©f~1  =  EUT. 

The  linear  dynamic  texture  model  is  found  by  minimizing  the  objective 

min^,  ©2  —  1  2F. 


Note  that 


where 


®2  »  [©"_1>0/v]  =  =  ©f"1  [Z)2l  VZ~1Utxn]  , 


D2  = 


0 

In 


-J 


From  (1.13)  we  have  VE =  s,  and  combining  this  with  (1.33, 1.51)  gives 

©^  =  ef-'s 


(1.31) 

(1.32) 

(1.33) 

(1.34) 

(1.35) 


The  optimal  solution  of  (1.32)  is  given  in  closed  form  by 

A  =  0^©f-lT  ’©f-1©jv-lT] _1  (1.36) 

=  e^'se?-'7,  [©"-1©"-17]'1 

=  G^-'SV E  [EV^VE]-1 

=  ©f-15VE_1,  (1.37) 


which  is  the  same  as  (1.25)  showing  that  the  dynamic  texture  model  is  the  same  as  the  linear  model 
used  in  DMD.  By  application  of  the  SVD  (1.10),  another  way  of  writing  (1.36)  is  simply 


A  =  ©^ [EVrVE] 

=  ©£'©(v-iTa_1  =  s. 


-i 
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(1.38) 

(1.39) 


Following  [23],  the  sample  input  noise  covariance  Q  is  estimated  by 


where  t )»  =  6l+\  —  Adi.  The  driving  noise  input  matrix  B  is  estimated  such  that 

BBt  =  Q. 


(1.40) 


(1.41) 


1.1.2  Nonlinear  Model  Identification  from  Spectral  Analysis  of  Markov  Operator 

We  provide  a  numerical  approach  to  estimating  nonlinear  stochastic  dynamic  models  from  time-series 
data.  After  possible  dimensional  reduction,  time-series  data  can  be  used  to  construct  an  empirical 
Markov  model.  Spectral  analysis  of  the  Markov  model  is  then  carried  out  to  detect  the  presence  of 
complex  limit  cycling,  almost  invariant,  and  bi-stable  behavior  in  the  model.  Model  parameters  are 
expressed  as  a  linear  combination  of  basis  functions  over  the  phase  space.  A  least  squares  minimization 
is  used  to  fit  the  basis  function  coefficients  in  order  to  match  the  spectral  properties  of  the  respective 
Markov  operators.  The  approach  is  demonstrated  on  the  estimation  of  a  nonlinear  stochastic  model 
describing  combustion  oscillation  data. 

We  provide  an  approach  for  estimating  nonlinear  dynamic  models,  possibly  driven  by  noise.  The 
estimation  approach  is  based  on  comparing  the  spectral  properties  of  the  empirically  constructed 
Markov  operator  with  the  model-based  Markov  operator.  A  nonlinear  model  is  fit  such  that  its  as¬ 
sociated  Markov  operator  has  similar  spectral  properties  as  the  empirical  Terms  in  the  stochastic 
differential  equation  model  are  estimated  by  a  linear  combination  of  basis  functions.  The  coefficients 
appear  in  the  numerical  approximation  of  the  Markov  model,  and  are  fit  using  least  squares  mini¬ 
mization.  Model  validation  is  motivated  spectral  methods  developed  in  [53]  [51]  for  the  comparison  of 

dynamical  systems. 

Spectral  method  for  analysis  and  comparison 

In  this  section,  we  describe  spectral  methods  for  the  analysis  and  comparison  of  the  dynamical  systems. 
The  material  for  this  section  is  taken  from  [20]  [53]  [51].  Also,  see  [45]  for  an  introduction  to  these 
concepts. 

Consider  the  stochastic  dynamical  system 

—  =  b(x)  +  <t{x)S,  xeXeRd,  (1.42) 

or  its  discrete  time  equivalent, 

xk+i  =  T(aj^.^fc)  (1.43) 

where  each  Xk  €  X  C  is  the  state  vector  and  €  U  is  sequence  of  i.i.d.  random  noise.  Associated 
with  T  is  a  stochastic  transition  function  p(x.  A),  which  gives  the  transition  probability  to  jump  from 
point  x  €  X  to  set  A  €  B(X),  where  B{X)  is  the  Borel  sigma  algebra  of  X.  For  deterministic  dynamics 
i.e.,  when  =  0,  we  have  p(x,  A)  =  ST^(A),  where  6  is  the  Dirac  delta  measure.  Stochastic  transition 


function  can  be  used  to  define  two  linear  transfer  operators  called  as  Perron- Frobenius  and  Koopman 
operators.  Here  we  consider  the  finite  dimensional  approximation  of  the  P-F  operator.  To  do  this  we 
consider  the  finite  partition  of  the  state  space  X  i.e., 


X  =  {Du...,Dm) 


(1.44) 


such  that  Di  fl  Dj  =  0  for  i  ^  j  and  =  X.  P-F  operator  on  the  finite  dimensional  vector  space 

Rm  can  be  represented  by  a  matrix  P  :  Rm  — ►  Rm  as  follows: 


Pil  = 


Mz.(r-1(Dj)nA) 

fitiDi) 


i,j  =  1, m 


(1.45) 


The  resulting  matrix  is  non-negative  and  because  T  :  Dt  — *  X,  £7=i  Pa  =  i  i.e.,  P  is  a  Markov  or  a 
row-stochastic  matrix. 

Important  complex  dynamical  features  of  the  dynamical  system  T  can  be  captured  using  its  Markov 
matrix  P.  For  example  long  term  or  asymptotic  behavior  of  the  dynamical  system  T  is  captured  by 
the  invariant  measure  or  more  appropriately  physically  relevant  measure.  Finite  dimensional  approx¬ 
imation  of  the  invariant  measure  or  the  outer  approximation  to  the  support  of  the  invariant  measure 
can  be  obtain  from  the  left  eigenvector  of  Markov  matrix  P  with  eigenvalue  one.  i.e., 


(xP  =  1  •  fJL 


Similarly  the  presence  of  periodic  or  limit  cycling  behavior  in  T  can  be  captured  by  the  complex 
unitary  spectrum  and  the  corresponding  eigenvectors  of  P.  Moreover  if  the  Markov  matrix  P  has  real 
eigenvalue  close  to  one  then  it  is  the  indicator  for  the  presence  of  almost  invariant  or  bistable  behavior 
in  the  dynamical  system  T.  For  more  detail  on  this  topic  refer  to  [19]  [20]. 

We  outline  an  approach  to  numerically  approximate  the  stochastic  dynamical  system  (1.42)  based 
on  the  empirically  obtained  Markov  matrix  P. 

Set 

d 

aij(x)  =  '52<7ik(x)crjk(x)  (1.46) 

k=  1 

Under  certain  regularity  conditions  [45],  the  evolution  of  the  density,  p  under  (1.42)  satisfies  the 
Fokker-Planck  equation, 


dp  __  1  v (i2  ,  .  d  ( 

dt  ~  2  dxidx/aiiP^  '  dxi  iP^ 

t  J=1  J  t=l 

•=Tp,  t>  0,x€Kd.  (1-47) 

Consider  the  finite-dimensional  (discrete-space)  approximation  of  T  by  discretizing  the  underlying 
phase  space  X  with  the  partition  X  following  (1.44).  Therefore,  the  distribution  p  is  approximated 
by  a  finiite-dimensional  vector  u(x)  €  Rm  and  let 

U  =  diag{u}  G  Rmxm,  (1.48) 
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where  diag{u}  is  a  diagonal  matrix  whose  entries  are  the  elements  of  the  vector  u.  Similarly,  the 
operators  are  approximated  by  their  finite-dimensional  matrices: 


fsf,  F  €  Rmxm. 

Continuing  this  way,  write  the  discretized  approximations  of  Oij'(x)  and  fei(x),  respectively, 

°y(x)  «  Aij(x)  e  Rm,  bi(x)  w  Bi(x)  6  Rm. 

Similarly  we  write  the  discretized  differential  operators 


1  9 2  «  Dfj  e  Rmxm,  •£—  «  £>i  €  R” 

OXi 


2  dxidij 


Next,  define  the  matrices: 


D2(U) 

D\U) 

D(U) 


=  [D2nU,...,D?jU,...,D2ddU]  GRm 
=  [DiU, . . . ,  DiU, . . . ,  DdU\  €  Rmxmd 
=  [D2(U),Dl(U)]  GR rnxim^+md) 


xmd2 


and 


‘  An  ‘ 

'  Bx  ‘ 

A  := 

Aij 

G  B:= 

Bj 

.  Add  . 

.  Bd  _ 

J )TTld 


C  := 


>1 

B 


j^mdP+rruf) 


(1.49) 


(1.50) 


(1.51) 

(1.52) 

(1.53) 


(1.54) 


The  infinite-dimensional  Fokker-Planck  operator  is  approximated  in  finite  dimensions  by 


Tp  as  Fu  =  D(U)C.  (1.55) 

We  next  write  (1.47)  in  discrete  time,  with 

du  ^  uj+i  —  ut 

dt  ~  It  ’ 

where  6t  is  the  time  step  which  the  time-series  data  was  obtained.  Substituting  this  into  the  left  side 
of  (1.47)  results  in  the  Markov  matrix  appearing  equation  (1.45), 

ut+i  =  ut  +  6tFut 

=  [ut  +  6tD(U)C] 

:=  Put.  (1.56) 
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We  empirically  obtain  P  directly  from  data  by  discretizing  the  underlying  phase  space  (Rrf)  to 
obtain  this  discrete-space  approximation.  The  spectrum  of  P  is  given  by  PV  =  VA,  where 


V  =  [v0,  t>i,  V2,  ■  ■ . ,  vm] ,  A  =  diag( A0,  Ai, . . . ,  Am). 

Under  suitable  conditions  there  exists  a  steady  distribution  t»o  where  Ao  =  1. 
Consider  a  single  eigenvector  v  £  Rm  with  eigenvalue  A.  We  have 


and  after  rearranging, 


Pv  =  [v  +  6tD(V)C]  =  Av, 


D(V)C  =  :=  w. 

dt 


(1.57) 


(1.58) 


(1.59) 


This  can  be  solved  for  C  which  will  give  the  functional  forms  of  6(x)  and  a(x )  appearing  in  (1.42).  Of 
course,  for  a  single  eigenvector,  there  are  multiple  solutions.  We  solve  (1.59)  for  multiple  eigenpairs; 
the  first  k  for  example.  We  have 


(1.60) 


Again,  for  k  small,  there  axe  multiple  solutions.  We  next  restrict  C  to  be  a  linear  combination  of  k 
basis  functions.  Let  {<pj}  be  a  set  of  basis  functions  where  <j>j  :  — >  R,  j  —  0, 1, 2, . . . ,  n  —  1,  such  as 
Hermite  polynomials.  Define  the  matrix  consisting  of  these  basis  functions  as  column  vectors: 


D(Vo) 

Wo 

c  = 

.  D{vk. i)  . 

.  wk-l  . 

We  rewrite  (1.54)  in  terms  of  these  basis  functions.  For  each  ifj  take 


and 


q  := 


Aij 

=  $Oij, 

a{j  €  Rn, 

Bi 

^5 

& 

II 

ft  6  Rn, 

an 

'ft  ' 

€Km|3) 

0:= 

ft 

.  add  . 

1 _ 

(1.61) 

(1.62) 

(1.63) 


e  R™*, 


c  := 


a 

0 


j^(nd2+Tuf) 


As  in  (1.51-1.53),  define  the  differential  operators  restricted  to  these  basis  functions: 
D%{U)  :=  [D2uU$,  DljUQ, . . . ,  DhUQ]  €  RmXnd2 

:=  [DiU$, . . . ,  DiU$ . DdU$]  €  Rmxnd 

:=  [D%(U),D\,(U)\  e  RmxM2+n<0. 


Di(U) 

D*{U) 


(1.64) 

(1.65) 

(1.66) 
(1.67) 
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Rewrite  equation  (1.60)  as 


:= 


DM) 


c  = 


wq 


Wk- 1 


:=Wk, 


(1.68) 


DM- 1)  . 

where  £>*,*  6  c  6  R("<**+nd) ,  Wk  €  Rfcm. 

We  solve  for  c  through  least  squares: 

min  [D<pikc  -  Wk]*  [D^kc  -  Wfc] ,  (1.69) 

C 

which  results  in 

•  c=[DlkD<s>tk]-1DitkWk-  (1-7°) 

We  would  like  to  compare  the  eigenfunctions  of  the  approximated  model  with  the  eigenfunctions 
appearing  in  equation  (1.57).  We  compute  the  approximate  Markov  matrix: 

d  d 

D(c)  :=  £  D% diag{$ati}  +  £  D^diag{<W?t}.  (1-71) 

ij=l  i=l 

The  approximate  Markov  matrix  is  then 

P  =  I  +  6tD{c).  (1.72) 

We  then  compare  the  spectrum  of  D{c)  with  the  spectrum  in  (1.57). 


Application  Example 

In  this  section,  we  construct  a  reduced  order  model  describing  nonlinear  oscillations  of  flame  dynamics. 
High-speed  video  data  was  obtained  from  the  UTRC  combustion  rig  described  in  [73],  and  proper 
orthogonal  decomposition  (POD)  was  used  for  dimension-reduction.  We  briefly  describe  how  the 
Markov  matrix  was  empirically  constructed. 


Data  Reduction  using  POD  Modes 


Prom  the  image  data,  the  mean  field  was  removed  from  each  of  the  images  and  the  POD  modes  were 
computed.  It  was  found  that  the  first  two  POD  modes  account  for  more  than  80%  of  the  energy  found 
in  the  data  set.  Hence,  we  consider  a  two  dimensional  state  space  x  €  R2.  The  two-dimensional  time 
series  is  constructed  via 


Xk  = 


Xl,k 


(<h,yk)  J 


(1.73) 


where  yk  denotes  the  A>th  image  with  the  mean  field  subtracted. 

The  time  series  of  POD  coefficients,  given  by  {x*},  were  obtained  from  projecting  the  original  data 
onto  the  POD  modes.  The  resulting  phase  space  (plotting  x\^  vs.  X2,jt,  for  k  =  1, . . .  ,iV)  is  shown 
in  figure  (1.1).  The  phase  portrait  shows  a  noisy  limit  cycle  where  the  density  of  points  is  clearly 
non-uniform.  The  rotational  speed  of  the  limit-cycle  varies  depending  on  the  state.  This  indicates 
that  a  nonlinear  model  is  necessary  to  match  this  time  series. 
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Projected  Phase  Space 


Figure  1.1:  The  phase  space  of  the  coefficients  resulting  from  projection  of  the  data  onto  the  first  two 
POD  modes 

Empirical  Markov  Model 

A  Markov  model  was  computed  from  this  time  series  data  through  equation  (1.45).  The  result  is  a 
non-reversible  Markov  model  with  many  eigenvalues,  as  shown  in  Figure  1.2.  Note  that  the  eigenvalues 
can  be  collapsed  toward  the  origin  by  taking  powers  of  P.  The  eigenvector  corresponding  to  the  unit 
eigenvalue,  shown  in  Figure  1.3  confirms  the  steady  distribution  of  the  trajectories. 


Figure  1.2:  The  eigenvalues  of  the  Markov  model. 

The  phase  of  the  eigenvector  associated  with  the  2nd  eigenvalue  is  shown  in  Figure  1.4.  It  reveals 
the  oscillatory  nature  of  the  dynamics. 

Nonlinear  Model  Extraction 

In  this  section  we  develop  a  second  order  stochastic  differential  equation  model  of  the  form 

ii  =  bi(xi,x2)  +  tfnfi 

±2  =  ^2(^1,  X2)  +  <722^2  (1*74) 

to  capture  the  essential  dynamical  behavior  of  the  reduced  set  of  data  from  the  previous  section. 
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magnitude  of  first  eigenvector 


Figure  1.3:  The  first  eigenvector  shows  the  invariant  distribution 


mi 


Figure  1.4:  Magnitude(left)  and  phase(right)  of  the  second  eigenvector 
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Before  going  into  the  details  of  the  model  we  would  like  to  summarize  some  of  the  key  dynamical 
features  of  the  reduced  set  of  data  as  captured  by  phase  portrait  in  Figure  1.1. 

•  Dynamics  consist  of  stable  limit  cycle  where  the  motion  along  the  limit  cycle  is  in  clockwise 
direction. 

•  The  limit  cycle  is  parameterized  by  angle  9.  Speed  along  the  limit  cycle  is  nonuniform,  speed  of 
the  limit  cycle  is  less  for  9  G  [0,  —  \ ]  compared  to  other  value  of  6 . 

We  would  like  the  model  to  capture  this  essential  dynamical  behavior  along  with  the  amplitude  of 
the  limit  cycle  and  the  average  speed  or  the  frequency  of  the  limit  cycle.  We  choose  a  reduced  set  of 
eigenvalues  closest  to  the  unit  circle  to  approximate  the  Markov  model.  Due  to  the  cyclic  nature  of 
the  data,  we  use  basis  functions  $  in  radial  coordinates,  expressed  as  separable  functions  in  r  and  6 

<£*,0(x1>x2)  =  r>C 
<i>k,2j{xi,x2)  =  rfccos  (j9) 

<Pk,2j+i(xi,x2)  =  rfcsin(j0), 

*  =  0,1,2,...,  i  =  0,1,2,:.; 

where 

x\  =  r  cos  9 
X2  =  r  sin  9 

The  least  squares  fit  (1.70)  resulted  in  basis  function  coefficients  producing  the  terms  appearing 
on  the  right  hand  side  of  (1.74).  The  functions  &i(xi,£2)  and  &2(xi, X2)  appear  in  Figures  1.5  and  1.6, 
respectively.  The  resulting  functions  crn{x  1,0:2)  and  <722(2:1, X2)  are  qualitatively  similar. 


-150  -100  -50  0  50 


Figure  1.5:  The  estimated  61(2:1,  £2)  appearing  in  (1.74). 

The  resulting  approximate  eigenfunctions  of  the  estimated  Markov  matrix  are  shown  in  Figure  1.7 
and  1.8.  Figure  1.7  shows  the  approximate  invariant  distribution  which  closely  resembles  that  shown 
in  Figure  1.3.  Similarly,  the  second  complex-valued  approximate  eigenfunction  is  shown  in  Figure  1.8 
which  closely  resembles  the  second  eigenfunction  shown  in  Figure  1.4.  The  match  in  the  eigenfunction 
indicates  a  good  match  between  the  model  and  the  data  in  terms  of  long  term  dynamics. 
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Figure  1.6:  The  estimated  62(^1, £2)  appearing  in  (1.74). 


Figure  1.7:  The  approximate  first  eigenvector  of  the  estimated  Markov  matrix 
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Figure  1.8:  The  approximate  magnitude(left)  and  phase(right)  of  the  second  eigenvector  of  the  esti¬ 
mated  Markov  matrix 
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1.1.3  Analysis  of  Complex  Spectra  and  Metastability 

The  purpose  of  this  paper  is  to  develop  methods  for  model  reduction  for  diffusion  processes  that 
exhibit  cyclic  behavior.  For  this  purpose  we  extend  techniques  based  on  the  spectral  theory  of  Markov 
processes  to  the  case  of  complex  spectra.  The  main  idea  is  to  augment  the  state  process  for  the 
diffusion  with  a  clock  process.  For  each  complex  eigenvalue  for  the  original  diffusion  there  exists  a 
real  eigenvalue  for  the  augmented  process.  Results  concerning  metastability  (or  quasi-stationarity)  axe 
then  applied  to  the  augmented  process.  For  the  special  case  of  a  linear  diffusion  in  two  dimensions, 
this  is  analogous  to  analyzing  the  process  in  a  rotating  coordinate  frame.  The  results  are  illustrated 
through  a  linear  diffusion,  and  an  empirical  model  of  combustion  dynamics. 

Extensions  of  the  classical  Wentzell-Freidlin  theory  for  model  reduction  have  appeared  in  numerous 
papers  over  the  past  decade.  Much  of  this  work  has  concerned  Markov  processes  that  are  reversible 
[59,  21,  13,  37,  14,  15].  The  goal  in  these  papers  is  to  understand  the  statistics  of  exit  times  from  a 
given  subset  of  the  state  space. 

Some  results  for  non-reversible  Markov  chains  are  available.  Fill’s  paper  [29]  extends  the  convergence- 
rate  bound  of  Diaconnis  and  Stroock  [22]  to  non-reversible  Markov  chains.  For  this  purpose  the 
transition  matrix  is  replaced  by  its  symmetrization,  and  the  rate  of  convergence  is  bounded  by  the 
eigenvalues  of  the  resulting  self-adjoint  matrix.  These  ideas  are  the  basis  of  [38]  that  establishes  exit 
time  statistics  from  a  set  for  a  discrete- time  non-reversible  Markov  chain. 

Extensions  of  Wentzell-Freidlin  theory  to  non-reversible  processes  appeared  for  the  first  time  in 
[38].  The  foundation  of  this  paper  is  the  theory  of  quasi-stationarity,  building  on  the  work  of  [27]. 
The  main  idea  of  [38]  can  be  summarized  as  follows:  Suppose  that  X  =  {X(£)  :  t  E  T}  is  a  diffusion 
process  evolving  on  X  =  with  transition  semigroup  denoted  {Pl  :  t  G  T}.  We  say  that  A  is  an 
eigenvalue  with  (non-zero)  eigenfunction  h  if  for  each 

Plh  =  eMh 

Suppose  that  A  is  real  and  negative.  In  this  case  we  can  assume  that  h  is  also  real-valued,  and  we  also 
assume  that  it  is  continuous.  We  would  like  to  consider  Doob’s  h-transform,  Pl  where 

Ig  is  the  multiplication  operator:  For  each  x  £  X  and  A  C  X  we  have  Ig(x,A)  =  p(x)I{x  €  A}.  The 
A-transform,  like  importance-sampling,  is  intended  to  lead  to  a  new  Markov  model  whose  properties 
provide  insight  into  the  problem  of  interest.  Unfortunately  {P*}  is  not  a  valid  Markov  semigroup  since 
A  may  take  on  negative  values.  Instead  we  consider  the  following  restricted  definition. 

Let  M  denote  a  connected  component  of  the  set  {x  :  h(x)  >  0}.  We  let  T#  =  inf(£  >  0  :  Y(t)  €  Mc), 
and  for  t  6  T  denote  =  t  A  T..  The  twisted  semi-group  is  defined  for  each  t  €  1\  x  €  X,  and  A  G  B 
(i.e.  A  Borel  measurable)  via, 

Pf(x,  A)  :=  [e-At-/i(X(t.))I{*(<.)  €  A}]  (1.75) 

Under  general  conditions,  it  is  shown  in  [38]  that  the  twisted  semi-group  corresponds  to  a  diffusion 
process  on  M  that  is  exponentially  ergodic.  Exponential  ergodicity  of  the  twisted  process  then  implies 
a  form  of  quasi-stationarity,  and  from  this  it  follows  that  the  exit  time  from  M  is  approximately 
exponentially  distributed  with  parameter  |A|. 
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The  inspiration  for  consideration  of  the  twisted  process  was  the  work  of  [27],  and  techniques  from 
the  large  deviations  analysis  contained  in  [5,  41]. 

The  main  result  of  this  paper  is  the  extension  of  the  results  of  [38]  to  the  case  in  which  A  €  C  is 
complex.  The  main  idea  is  to  augment  the  state  process  for  the  diffusion  with  a  clock  process.  For  each 
complex  eigenvalue  for  the  original  diffusion  there  exists  a  real  eigenvalue  for  the  augmented  process. 
Results  concerning  metastability  contained  in  [38]  are  then  applied  to  the  augmented  process. 

It  is  assumed  that  X  =  (AT(t)  :  t  G  T}  is  a  diffusion  process  evolving  on  X  =  Rd,  with  transition 
semigroup  denoted  {Pl  :  t  G  T}.  Letting  u  denote  the  drift,  and  E  the  covariance  matrix,  the 
differential  generator  (see  e.g.  [45])  is  defined  for  C1  functions  h:  X  — *  C  by  Vh  (x):= 

E-We‘( +  (i.re) 

i  ij  J 


or,  in  more  compact  notation, 

V  =  u  •  V  +  ^trace  (EA) 

For  each  (3  >  0  the  resolvent  kernel  is  given  as  the  Laplace  transform, 

R0:=  [°°  e-ftp*  dt. 

Jo 


(1.77) 


We  write  R  :=  Rp  when  /?  =  1. 

It  is  assumed  as  in  [38]  that  the  diffusion  is  V -uniformly  ergodic:  For  a  probability  measure  v  on 
B ,  some  constants  b  <  oo  and  F  >  0,  a  function  s  :  X  — *  [0,  oo),  and  a  V  :  X  — >  [1,  oo): 


VV  <  ~rv  +  bs 
R  >  s  <g>  v. 


(V4) 


The  second  inequality  in  (V4)  means  that  the  function  s  and  probability  measure  v  are  small  This 
terminology  and  the  outer  product  notation  axe  taken  from  [57].  This  ‘smallness  assumption5  is 
equivalently  expressed. 

R(x,A)  >  s(x)i/(j4),  x  e  X,  A  €  B. 

Suppose  that  V  has  a  complex  eigenvalue  A,  which  we  write  as 


A  =  — T  +  id 

with  T  >  0,  and  d  ^  0,  with  associated  eigenvector  h.  Consider  the  clock  process  defined  by, 

$(t)  =  $(0)e^,  t  >  0,  (1.78) 

with  initial  condition  restricted  to  the  unit  circle  in  C,  which  is  denoted  U.  The  clock  process  is 
Markov,  as  is  the  bivariate  process, 


Y(t)  = 


t  >  0. 
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In  fact  Y  is  a  diffusion  on  Y  =  X  x  U  whose  covariance  matrix  for  y  is  given  by, 


Er(j/):=diag(E(x),0).  (1.79) 

Throughout  the  paper  we  adopt  the  notation  y  =  ( x ,  4>)  for  y  €  Y.  with  x  €  X,  <j>  G  U. 

We  define  for  each  real  /?  G  R  the  real- valued  function, 

90(y)  =  Re  ((el0/<f>)h{x)),  y  =  (x,  <j>)  G  Y.  (1.80) 

Proposition  1.1.1  For  each  ft  G  R  the  function  gp  is  an  eigenfunction  for  the  process  Y ,  with 
eigenvalue  A y  =  — I\ 

Proof.  The  differential  generator  for  X  can  be  extended  in  the  obvious  way  to  Y.  Given  the 
simple  dynamics  of  we  have  for  any  function  /:  U  — ►  C, 

Vf(<f>)  =  iW'(4> ) 

With  f(</>)  =  I /<p  the  eigenfunction  equation  holds, 

vf  (*)  =  =  -i^m,  <f>  g  u. 

Hence  the  generator  applied  to  gp  gives, 

Vg0(y)  =  Re((ei0/<f>)Vh(x)) 

+  Re  (x)) 

=  Re  ((eif}/<t>)Ah{x) 

4 — / <}>)h  (x)) 

=  -r g0(y),  y€  Y. 


The  twisted  process 

To  define  the  twisted  process  we  fix  (3  =  0  in  the  definition  (1.80),  and  let  M  denote  a  connected 
component  of  {y  :  po(y)  >0}.  It  is  assumed  that  this  set  has  nice  topological  properties:  M  is  equal 
to  the  closure  of  its  interior.  Following  [38],  we  define  T#  =  inf(t  >  0  :  Y(t)  G  Mc),  and  the  associated 
twisted  process  as  follows: 

The  twisted  process  is  the  Markov  process  Y  with  state  space  M  whose  semigroup  is  defined 
using  (1.75)  based  on  the  eigenfunction  g$.  Equivalently,  for  each  /  G  Lco(M),  and  any  x  G  M, 
Psf(y):=Ey[f(Y(s))]  = 

^E*[p0(Y(s  A  T.))f(Y(s  A  T.))  exp((s  A  T.)T)] 

The  twisted  process  has  a  generator  defined  for  C 2  functions  / :  Y  — »  C  by, 

Vf  =  g^Vigof)  +  Tf.  (1.81) 
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Two  key  assumptions  are  imposed  in  [38]:  First,  that  the  diffusion  is  hypoelliptic  (which  is  used 
to  conclude  that  the  resolvent  possesses  a  density  with  respect  to  Lebesgue  measure).  Second,  it  is 
assumed  that  the  gradient  of  the  eigenfunction  does  not  vanish  on  the  boundary  of  M.  The  gradient 
assumption  is  maintained  here.  To  ensure  that  Y  is  hypoelliptic  we  assume  that  X  is  elliptic ,  meaning 
that  its  covariance  is  strictly  positive.  These  assumptions  are  collected  together  as  follows: 

E(x)  >  0, 

V,so(y)  =  Re(< r'Vh(x))  ?  0 ,Vy  €  9M.  (1.82) 

It  is  not  hard  to  see  that  the  assumption  (1.82)  always  fails  when  X  is  a  diffusion  in  one-dimension. 
We  see  in  the  next  section  that  it  does  hold  in  many  examples,  such  as  the  linear  diffusion  in  two  or 
more  dimensions. 

The  following  result  is  a  consequence  of  Theorem  3.7  of  [38].  The  reader  is  referred  to  this  paper 
for  a  precise  definition  of  metastability  —  Its  main  conceptual  conclusion  is  that  the  exit  time  T.  is 
approximately  e-xponentially  distributed,  and  that  the  process  ‘almost’  reaches  a  ‘local’  steady-state 
prior  to  exiting  M. 


Theorem  1.1.2  Assume  that  (V4)  is  also  satisfied  for  a  continuous  function  V :  X  — >  [l,oo).  Suppose 
that  h  is  an  eigenfunction  with  complex  eigenvalue  A  =  —  T-f  id  satisfying  the  following  conditions: 

(a)  0  <  T  <  r. 

(b)  go(y)  >  0  for  all  x  e  M,  and  go(x)  =  0  for  x  €  dM  ;=  M  \  M. 

(c)  Condition  (1.82)  holds .  Consequently ,  fory  €  dM ,  (V5o(^))T2y(y)(Vffo(y))  >  0. 

(d)  Kn  :=  {x  £  X  :  V(x)  <  is  a  compact  subset  ofX  for  each  n  >  1. 

Then , 

(i)  The  escape-time  from  M  for  the  twisted  process  is  infinite  a.s .  for  y(0)  =  y  €  M; 

(ii)  The  twisted  process  is  V\  -uniformly  ergodic  with  V\  (y)  =  V(x)/go{y),  y  6  Y. 


(iii)  The  set  M  is  both  metastable  and  V -metastable,  with  exit  rate  T(M)  =  IV(M)  =■  T.  In 
particular, 


E[eer*] 


ife>  T 
otherwise. 


□ 


The  proof  of  Theorem  1.1.2  amounts  to  establishing  a  version  of  (V4)  for  the  twisted  process.  We 
can  follow  the  same  steps  as  in  [38]  to  construct  the  required  Lyapunov  function. 

For  a  given  0  <  a  <  1  write 


Vi  :=  So  1 V,  V2  :=  So  .  and  V  :=  V,  +  V2  . 
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We  denote  Go  =  \og(go),  where  go  is  the  eigenfunction  for  Y.  From  (V4)  and  the  eigenvector  equation 
we  have, 

VV,  =  \I~'VIgo+TI}g^V 
=  g^lVV  +  TV] 

<  _(r  -  r)Vj  +  bgols 

t>V2  =  [Ig^VIgo  +  ri]go1+a 
=  9o  1[®5o  +  ^5ol 
=  Q^_1[r  -1(1-  a)VG0rSyVGo]. 

Following  arguments  in  [38],  we  obtain  a  version  of  (V4)  for  the  twisted  process:  For  a  finite  constant 
&o,  and  a  compact  set  S  C  M, 

w^-itr-rjy+Ms. 

Examples 

We  discuss  an  analytic  example  as  well  as  an  example  motivated  by  an  empirical  model  of  limit-cycling 
combustion  dynamics. 

Ornstein-Uhlenbeck  process 
Consider  the  Ornstein-Uhlenbeck  process, 

dX  ( t )  =  AX{t)  dt  +  dW{t),  (1.83) 

where  W  is  a  full-rank  Gaussian  process.  Suppose  that  A  is  a  complex  eigenvalue,  and  v  a  (non- zero) 
left-eigenvector  for  A ,  satisfying 

ATv  —  Au. 

The  generator  for  X  shares  this  eigenvalue,  and  the  function  h{x)  =  vrx  is  an  eigenfunction: 

Vh(x)  =  (Ax)TX7h(x)  +  ^trace(EA/i(x)) 

=  xTATv  =  A  h(x). 

We  now  check  to  see  if  (1.82)  is  satisfied.  We  have, 

Vi9o(y)  =  Re  y  =  ( X ,  (j))  €  Y. 

This  is  zero  if  and  only  if  Re  (0~1ua:)  =  0  for  each  k  =  1, . . . ,  n.  If  this  holds  for  some  0  6  U,  it  then 
follows  that  v *  =  i<p~lv  is  a  purely  real  eigenvector  for  A,  which  is  impossible  since  A  is  complex.  We 
conclude  that  (1.82)  is  satisfied. 

Consider  the  two-dimensional  model  with 


22 


where  a  >  0.  The  matrix  A  possesses  a  pair  of  complex  eigenvalues  in  the  left-hand  complex  plane, 
satisfying  T  =  a: 

eig  (A)  =  —  a  ±  i. 

A  left  eigenvector  for  A  is  given  by  vT  =  [—1,  *],  which  gives 

Re(e-JV*(t))  =  cos(t)Xi(0  +  sin(£)X2(£). 


If  X (0)  satisfies  Re  (vTX (0))  >  0,  we  can  expect  that  Re(e-JtvrX(f))  >  0  for  a  period  of  time 
approximately  exponentially  distributed,  with  mean  1/a.  Applying  Theorem  1.1.2  we  conclude  that 
the  first  exit  time  T.  =  inf(£  >  0  :  Re  (e~:>lvTX(t))  =  0)  shares  the  following  property  with  the 
exponential  distribution: 

E[<^](=0°  if£2‘‘ 

I  <  oo  otherwise. 


□ 


Empirical  Model  of  Limit-Cycling  Combustion  Dynamics 

We  apply  the  analysis  to  a  Markov  model  describing  the  nonlinear  dynamics  of  limit-cycling  combustion 
oscillations.  The  data  was  obtained  from  an  experimental  combustion  rig  described  in  [73].  The  two- 
dimensional  phase  space  was  obtained  as  in  [?]  as  follows.  A  POD  analysis  was  done  on  the  temporal 
flame  images  and  the  data  wras  projected  on  to  the  first  two  dominant  POD  modes.  The  dynamics  of 
the  flame  data  projected  on  to  this  two-dimensional  space  is  shown  in  Figure  1.1.  The  phase  portrait 
shows  a  noisy  limit-cycle  where  the  direction  of  oscillation  is  in  the  clockwise  direction. 

A  discrete  time  Markov  model  was  constructed  for  the  dynamics  on  this  two-dimensional  space. 
The  eigenvalues  are  shown  in  Figure  1.9.  The  complex  eigenvalues  suggest  cyclic  behavior  and  a 
metastability  analysis  can  be  done  using  the  corresponding  eigenfunctions  as  described  in  the  previous 
sections. 
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Figure  1.9:  Eigenvalues  of  the  Markov  matrix  associated  with  the  combustion  dynamics  data  shown 
in  Figure  ??. 

We  describe  the  metastable  sets  associated  with  the  eigenvalues  shown  on  the  right  in  Figure 
1.9.  In  particular,  the  eigenvalue  at  A  :=  |A|e^  =  0.98  +  jfO.995  is  associated  with  an  eigenfunction 
that  varies  in  the  tangential  direction  and  has  no  radial  variation.  The  associated  eigenvector  h{x) 
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is  complex  as  shown  in  Figure  1.10.  We  take  the  clock  process  to  be  the  discrete  time  equivalent  of 
(1.78), 

<t>k  =  k  =  1, 2, . . . , 


where  ip  is  the  angle  of  the  eigenvalue  A.  Setting  </>$  =  1,  the  eigenfunction  of  the 
process  is 


g0(y)  =  Re  (e  lt/fkh{x))}  y  =  (x,  <p)  G  Y. 


associated  twisted 


Figure  1.10:  The  complex  eigenvector  h(x)  (magnitude  on  right,  phase  angle  on  left)  associated  with 
the  complex  eigenvalue  A  =  0.98  +  jO. 995  shown  in  Figure  1.9. 

The  plot  in  Figure  1.11  shows  the  sign  of  go{y)  for  different  phase-shifts  (i.e,  after  multiplication 
by  e~%^k  for  different  values  of  k.).  Note  how  the  sets  with  positive  support  and  negative  support 
rotate  around  the  phase  space  and  the  exit  time  marks  the  point  when  the  system  exits  one  of  these 
rotating  sets  (i.e.,  exhibits  a  phase-shift  in  its  oscillations). 


Figure  1.11:  The  sign  of  the  eigenfunction  with  a  complex  eigenvalue  close  to  the  unit  circle,  rotating 
with  incremental  phase-shifts  of  j  betwreen  0  and 

The  eigenvalue  A  =  0.89  is  purely  real  and  hence  has  a  purely  real  eigenvector  with  no  tangential 
variation,  but  variation  in  the  radial  direction.  The  sign  of  the  eigenvector  is  shown  in  Figure  1.12. 
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Since  the  eigenvalue  is  real,  this  eigenvector  is  not  associated  to  any  cyclic  behavior.  This  eigenvector 
and  the  related  exit  time  simply  indicates  when  the  system  moves  from  a  state  of  low  amplitude 
oscillation  to  high  amplitude  oscillation,  and  vice-versa. 


absevecS 


Figure  1.12:  The  sign  of  the  radial  eigenvector  of  the  Markov  matrix. 

Finally,  the  complex  eigenvalue  A  =  0.851H- jr0.99  has  an  eigenvector  exhibiting  both  tangential  and 
radial  components.  Note  again  how  the  metastable  set  rotates  around  the  phase-space,  as  indicated 
by  the  phase-shifted  sign  of  the  eigenvector  shown  in  Figure  1.13. 


Figure  1.13:  The  sign  of  the  eigenvector  with  tangential  and  radial  variation,  shown  rotating  with 
incremental  phase-shifts  of  f  between  0  and 

By  examination  of  the  magnitudes  of  the  eigenvalues,  the  eigenvectors  associated  with  these  three 
metastable  sets  have  decreasing  mean  exit  times.  This  is  intuitively  confirmed  by  the  fact  that  the 
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sets  become  increasingly  complicated.  A  hierarchy  of  such  sets  along  with  the  spectral  properties  of 
the  Markov  matrix  can  be  used  to  construct  a  reduced  order  model  of  the  measured  process  through 
techniques  described  in  [64]. 

We  have  presented  a  framework  for  analyzing  Markov  models  with  semi-rotational  dynamics  by 
considering  the  complex  spectra,  and  illustrated  the  approach  using  an  application  involving  limit¬ 
cycling  combustion  oscillations.  The  ultimate  goal  of  this  research  is  to  construct  low  order  models 
that  capture  essential  structure,  such  as  the  hidden  Markov  models  proposed  in  [38].  The  most 
interesting  open  problems  are  application  specific.  For  example,  can  we  justify  the  consideration  of 
a  two-dimensional  model  obtained  from  POD  coefficients?  If  not,  what  are  alternative  approaches  to 
treat  the  full-order  Maxkov  model? 

1.2  Model  Reduction 

1-2.1  Tangent  Space  Approach  to  Nonlinear  Model  Reduction  and  Identification 

In  this  section  we  discuss  a  novel  approach  for  model  reduction  of  nonlinear  systems  with  output 
measurements  based  on  the  analysis  of  linear  derivative  maps  [74].  With  every  nonlinear  system 
one  can  associate  a  linear  derivative  map  that  evolves  vectors  on  the  tangent  space.  At  each  point 
along  a  nominal  trajectory,  a  local  observability  gramian  defined  on  the  tangent  space  is  computed 
based  on  the  linear  perturbation  dynamics,  which  is  then  used  to  identify  a  balanced  local  coordinate 
system.  These  local  coordinates  can  be  patched  together  to  construct  the  global  coordinates  for  the 
reduced  order  representation  of  the  system.  A  computational  approach  is  described  for  the  empirical 
construction  of  the  gramian  on  the  tangent  space  and  for  the  alignment  of  the  local  coordinates  to 
obtain  the  global  coordinates.  Simulation  results  and  examples  are  presented  to  demonstrate  the 
application  of  the  proposed  method. 

Previous  methods  of  model  reduction  include,  for  example,  techniques  based  on  proper  orthogonal 
decomposition  (POD)  for  of  fluid  flows  problem  [11],  [62],  balanced  truncation  in  control  systems  [54], 
[65],  [43],  spectral  analysis  of  transfer  operators  [51],  [52],  and  fast  and  slow  manifold  decomposition 
[30]  in  dynamical  systems. 

Derivative  Gramian  and  Model  reduction 

Consider  the  model  reduction  problem  for  a  discrete  time  dynamical  system  with  output  measurement 
as  follows: 

Xjfc+1  =  f&k) 

Vk  =  h{xk)  (1.84) 

where  x  €  X  C  Rn  is  a  compact,  state  space  and  y  6  Y  C  Rm.  Associated  with  the  nonlinear  system 
is  a  linear  derivative  map  obtained  from  the  linearization  of  the  system  along  a  nominal  trajectory. 
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The  linear  derivative  map  evolves  vectors  on  the  tangent  space  and  is  defined  as  follows: 

*)k+l  =  =:  A{xk)Vk 

(j>k  =  Xk)Vk  —■  C[xk)rik  (1.85) 

where  A(xk)  :  TXlcX  -*  Tj^X  and  C(xk)  :  TXkX  -*  Th^Xk)Y.  Important  information  about  the  sys¬ 
tem  dynamics  can  be  obtained  from  the  linear  derivative  map.  For  example,  the  Lyapunov  exponent, 
which  can  be  thought  of  as  the  generalization  of  eigenvalues  from  linear  systems  to  nonlinear  systems, 
are  obtained  from  the  linear  derivative  map.  The  negative  (positive)  value  of  maximum  Lyapunov 
exponent  implies  exponential  convergence  (divergence)  of  nearby  system  trajectories.  Evaluation  of 
the  derivative  map  at  the  trivial  solution  is  widely  used  in  the  local  stability  analysis  of  the  trivial 
solution. 

Under  suitable  technical  conditions  [74],  for  a  given  point  x  on  the  nominal  trajectory  the  observ¬ 
ability  gramian 

oc  /  A: — 1  \  fc— 1 

Q(x)  =  E  n  C'(xk)C(xk)  n  A{Xj)  (1.86) 

fc=0  \j=0  /  j= 0 

is  well  defined.  It  is  easy  to  check  that  Q(x)  is  positive  semi-definite  and  hence  defines  a  pseudo¬ 
metric  on  the  tangent  space  at  x,  with  {£,  7i)q(x)  =  £Q(x)t?,  where  f  and  t\  are  vectors  belonging  to 

TxX.  Furthermore,  one  can  verify  that  Q(x()  =  X (xt)Q{xi+i)A{xt)  +C  (xi)C(xt)  where  x/+i  = 
f(xt),  which  resembles  equation  (7)  in  [42],  where  Lall  and  Beck  write  the  Lyapunov  inequality  for  a 
generalized  observability  gramian. 

Since  Q{x)  is  positive  semidefinite  we  know  that  there  exists  a  unitary  transformation  U (x)  as 
a  function  of  x  such  that  U  (x)U(x)  =  I  and  £(x)  =  U(x)Q(x)U' (x)  is  a  diagonal  matrix.  Now 
consider  the  coordinate  transformation  on  the  tangent  space  if  —  U(x)t],  t]  —  U'  (x)ip  so  that  xpk+i  = 
U{xk+\)A(xk)U' {xk)ij)k  :=  A(xk)ipk  and  <j>k  =  C{xk)U' (xk)vk-  The  derivative  gramian  in  the  new 
coordinates  can  be  written  as  Q(x/)  =  (xe~1)C(xk)'  C(xk)A(x^~l).  Now 

-4(x£-1)  =  A(xk-\)A(xk-2)  •  •  •  A(xi) 

=  U(xk)A(xk-i)A(xk-2)  ■  ■  ■  A{xe)U'(xt) 

=  UixJAixl-^U'ixt)  (1.87) 

Hence  we  have 

oo 

Q(xt)  =  [7(x/)E^,(^“1)C''(xfc)C(xfc)^(xJ-1)U'(x<) 
k=t 

=  U(xe)Q(xe)U\xi)  =  Z(xt) 

Note  that  the  local  coordinate  U[x)  defined  at  each  point  of  the  tangent  space  provides  a  diagonal 
decomposition  of  the  derivative  gramian  at  each  point  x.  One  would  like  to  patch  these  local  coordinate 
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to  construct  the  global  coordinate  on  the  phase  space.  Here  we  proceed  formally  and  assume  that 
such  global  coordinate  system  exists  i.e.,  z  =  T(x)  such  that  =  U(x).  In  the  next  subsection 
we  provide  a  computational  approach  for  aligning  these  local  coordinates  to  approximate  the  global 
coordinates. 

We  assume  that  in  the  diagonal  matrix  £(x)  has  the  following  structure. 


and  that  £1  >>  £2-  Let  z  =  T(x)  and  we  partition  the  state  z  based  on  the  partition  of  its  derivative 
gramian  £  i.e.,  z  =  ( zl,z 2)  =  (T\ (a:) , T2 (2$) )  =  T(x). 

It  turns  out  that  the  reduced  order  system  with  the  output  measurement 

4+i  =  TxU{T-\zl  0))) 

Vk  =  h(T~l{zlO))  (1.88) 

can  be  shown  to  have  non-positive  Lyapunov  exponents  with  diagonal  derivative  gramian  £i(£fc,0). 

The  computation  of  the  empirical  gramians  is  similar  in  spirit  to  [43]  and  [62]  but  with  a  key 
difference:  the  formulae  in  [43]  subtract  the  temporal  means  from  the  trajectories,  which  is  not 
applicable  in  the  current  setting.  Indeed,  for  the  empricial  gramian  based  on  the  derivative  mapping, 
we  must  track  the  trajectories  of  the  nominal  mean  states  and  outputs. 

Local  Tangent  Space  Alignment 

The  global  coordinate  system  T  computed  from  the  localized  coordinates  so  that  ^  U(x)  can  be 
approximated  through  an  alignment  procedure.  Denote  the  K -dimensional  vector  all  ones  by  e#.  We 
approximate  the  global  coordinates  T  using  the  initial  condition  data  points  in  {A/o  . . .  A/}v}  following 
the  procedure  of  [79].  Further  analysis  and  refinements  appear  in  [78,  76,  46].  Locally,  the  balancing 
coordinates  are  given  by  the  transformation  tyk  =  t/'(:rfc)A/i;.  Suppose  that  there  is  sufficient  overlap 
between  the  sets  of  initial  conditions  [78]  and  that  there  are  M  distinct  initial  condition  points. 
Furthermore,  suppose  that  the  partition  is  such  that  z\  6  In  other  wrords,  the  reduced  order 
model  will  be  of  dimension  d.  We  now  construct  the  global  coordinates  rk  for  k  =  1, . . . ,  N.  Write 
Tk  =  . . .  ,r*]  and  seek  to  minimize  Ek  =  Tk  ( I  —  —  Lk^k  for  each  k.  The  alignment 

matrices  Lk  can  be  determined  separately  by  use  of  the  Moore-Penrose  pseudo-inverse.  The  global 
coordinates  7).  are  approximated  by  minimizing 

E  lISfcllF  =  £  \\TkWk\\2F  =  HrSW’llJ.  =  IV  (TSWWtStTt)>  (1.89) 

k  k 

where  T  =  [rx, . . . ,  rM]  €  RdxM ,  S  =  [Si, . . . ,  SN]  €  where  Sk  €  RMxK  is  the  0-1  selector 

matrix  such  that  TSk  =  Tkt  and  W  =  diag  {W\, . . . ,  Wjy).  We  apply  the  constraint  that  TTt  =  /. 
The  rows  of  T  are  taken  to  be  the  eigenvectors  associated  with  the  next  smallest  2nd  through  d  +  1st 
eigenvalues.  The  reader  is  directed  to  [79]  for  further  details  on  this  computation. 
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Example:  Nonlinear  System  Reduction 

We  illustrate  the  computation  of  the  empirical  gramians  and  the  global  coordinates.  Consider  the 
controlled  nonlinear  system 

nt+i  =  rk  4-  0.1  (R  -  rfc)  (1  +  cos(lO0)) 

#fc+i  =  0k  +  A  4-  (uk  —  0k) 

Zk+i  =  zk  (1  +  O  O1(1+COS(0))) 

Ufc  =  Ufc_ i  +  A,  uo  =  0 

yk  =  [rjfc  cos (0k),  rk  sin(0fc),  zk] .  (1.90) 

This  system  has  equilibrium  rk  —  R,  6k  =  uk,  zk  =  0.  An  example  of  perturbation  and  nominal 
trajectories  with  A  =  ^  are  shown  in  figure  1.14.  Notice  that  the  ^-coordinate,  while  small  relative 
to  the  radius  of  R  =  1,  decays  slowly. 


Figure  1.14:  Nominal  trajectory  shown  by  red  circles  and  example  perturbation  trajectories  shown  by 
blue  dots  and  stars.  Initial  perturbations  are  in  the  +/  —  z  (dots)  and  +/  —  r  (stars)  directions. 

A  set  of  initial  conditions  shown  in  figure  1.15  (left)  was  used  to  compute  the  empirical  derivative 
gramians  and  the  global  balancing  coordinate  system.  Locally,  there  are  two  dominant  singular  values 
of  the  gramians,  based  on  the  simulation  results.  The  resulting  principle  component  vectors  (scaled  by 
the  singular  values  of  the  gramian)  associated  with  four  points  are  shown  in  figure  1.15  (right).  The 
vectors  show  that  there  are  two  dominant  component  directions  in  the  0,  z  plane.  Although,  geomet¬ 
rically  the  ^-coordinate  is  small  compared  to  the  0  coordinate  and  about  equal  to  the  r  component,  it 
is  dynamically  more  important  because  it  decays  slowly. 

The  principle  components  give  only  local  information  about  the  vector  field.  We  employ  the  align¬ 
ment  procedure  described  in  section  1.2.1  to  approximate  the  global  coordinates  T.  Although  we  have 
previously  discussed  the  local  principle  components  associated  with  the  derivative  dynamics  in  terms 
of  the  global  coordinate  system,  the  alignment  procedure  requires  no  a  priori  information  regarding 


(a)  Nominal  Trajectory 


(b)  Principal  Components 


Figure  1.15:  Left:  Nominal  trajectory  shown  by  red  circles  and  the  set  of  initial  conditions  used 
to  compute  the  empirical  gramians.  Right:  Local  principle  component  directions  computed  at  four 
different  points  on  the  nominal  trajectory. 


the  global  structure  of  the  principle  directions.  Figure  1.16  (left)  shows  the  resulting  approximation 
of  the  0  coordinate  based  on  the  global  alignment  of  local  principle  component  direction.  The  figure 
shows  (with  the  appropriate  phase  shift)  almost  a  1-1  relation  between  the  empirically  determined 
global  coordinate  and  the  actual  0  value  of  the  initial  condition. 

Figure  1.16  (right)  shows  a  similar  relation  between  the  empirically  determined  z  coordinate  and 
the  actual  z  values  of  the  initial  condition  points.  Again  (with  scaling  of  -1)  there  is  a  near  1-1 
correspondence.  This  example  illustrates  that  a  set  of  reduced  global  coordinates  can  be  approximated 
from  the  local  gramian  information,  with  no  a  priori  information  regarding  the  underlying  coordinate 
system  of  the  dynamic  model. 

Extensions  to  Dynamic  Texture  Models 

The  tangent  space  approximation  in  terms  of  local  PCA  and  the  computations  of  the  empirical  ob¬ 
servability  gramian  provides  a  framework  for  calibrating  hybrid  affine  or  nonlinear  dynamic  texture 
models,  which  is  the  subject  of  ongoing  work  at  UTRC  currently  funded  by  AFOSR.  Following  the 
methodology  for  identifying  balanced  stochastic  systems  in  [3]  we  can  use  the  tools  described  in  this 
section  to  identify  generalizations  to  the  dynamic  texture  models  originally  presented  in  [23].  Such 
generalizations  are  appropriate  for  video  data  that  exhibits  nonlinear  behavior  or  mode  switching 
which  can  be  modeled  by  a  piecewise  affine  [61]  dynamic  texture.  This  is  similar  in  spirit  to  the 
identification  methods  based  on  statistical  clustering  appearing  in  [28,  56]-  Based  on  the  local  PCA 
representation  of  the  data,  locally  affine  models  can  be  calibrated  based  on  standard  L*i  techniques 
[49].  The  model  switching  is  determined  by  relative  distances  between  the  model’s  current  state  and 
the  local  means  of  neighboring  data  patches.  The  concept  of  nonlinear  dynamic  textures  has  appeared 
in  [48,  4,  47,  77].  An  example  of  a  piecewise  affine  dynamic  texture  model  of  the  beach  video  sequence 
[4]  is  shown  in  Figure  1.17  (bottom  row),  where  it  is  compared  with  a  linear  dynamic  texture  model 
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(a)  6  coordinates  (b)  z  coordinates 


Figure  1.16:  Left:The  6  component  of  the  estimated  global  coordinates  T  compared  to  the  9  com¬ 
ponents  of  the  initial  condition  sets  Mq,  .  - .  ,A/jv-  Right:  The  z  component  of  the  estimated  global 
coordinates  T  compared  to  the  z  components  of  the  initial  condition  sets  JV/q,  •  •  • ,A/)v. 


(top  row).  The  top  row  shows  that  the  frames  of  the  linear  dynamic  texture  model  capture  the  bulk 
motions  of  the  waves  and  eventually  saturate.  The  piecewise  affine  model  frames  shown  in  the  bottom 
row  capture  the  different  dynamic  modes  of  the  waves  coming  inward  and  outward  at  the  same  time 
increasing  the  image  quality  and  avoiding  saturation  effects. 


1.2.2  Empirical  Controllability  Gramians  from  Trajectories 

In  this  section  we  provide  more  details  on  the  computation  of  empirical  gramians. 

Suppose  that  for  the  controlled  system  we  have  the  nominal  and  perturbation  trajectories  {zfc} 
and  {r)k},  respectively,  resulting  from  a  nominal  control  input  {u*}  and  an  impulse  input  v1.  We  can 
write  the  state  history  in  terms  of  the  impulse  input  at  k  =  0 


’  nh  ' 

i 

to 

H 

o 

'cT 

92,1 

^(xijR^o)^1 

nli 

A(x2)A(xi)B(xo)v1 

.  : 

.  •’ 

The  empirical  controllability  gramian  can  be  computed  from  a  collection  of  impulse  responses.  If 
we  apply  the  unit  impulse  u1  at  k  =  —j  then  the  response  is 


"  nzki " 

B(x-j)v 1 

n -l+2,\ 

A(x-j+i)B(x-j)vl 

n-U 3,1 

A(x-j+2)A(x-j+\)B(x-.j)vi 

• 

- 

(1.92) 
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(e)  Frame  20  (f)  Frame  40  (g)  Frame  60  (h)  FVame  80 

Figure  1.17:  Top  Row:  Linear  dynamic  texture  of  beach  sequence.  Bottom  Row:  Hybrid  affine 
dynamic  texture  of  beach  sequence. 


so  that 


'  V1J ' 

’  A(xo)  •  •  •  A(x-j)B(x-j)v1 

*72,1 

A(xi) . . .  A(x-.j)B{x-j)vl 

%,i 

A{x 2) . . .  A{x-j)B(x~j)vl 

Suppose  we  assemble  a  collection  of  impulse  inputs 

V:=[v\  ....  v*] 

and  their  corresponding  state  responses  if  the  impulse  inputs  were  applied  at  k  =  —j 


A(x0)  •  ■  -A(x-j)B(x-j) 
A{x{) . . .  A{x-j)B{x-j) 
A(x2)  ■  ■ .  A(x-j)B(x-j) 


V, 


(1.93) 


(1.94) 


(1.95) 


where  the  Af(x_j)mi*  =  is  the  state  response  at  time  m  with  respect  to  unit  impulse  input  vl 
applied  at  time  k  =  —j.  Denote  the  m-th  block  row  of  Af(x-j)  by  Af(x-j)m.  Assume  that  V  is 
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orthonormal. 


Y^-ih^'^-jh  =  B(xq)W'B'(xq)  +  A(x0)B(x.1)VVB,(x.l)A,(x0)  (1.96) 

5= 0 

+  A(xQ)A{x--l)B(x-2)VVB,{x-'l)A'{x-i)A\x0)  +  •  •  • 

=  B(xo)B'(x0 )  +  (1-97) 

+  A{xo)A{x-i)B(x-i)B'  (x_2M'  {x-i)A'  (xo)  +  . . . 

=  B(x0)B'(x0)  =  V{xx).  (1.98) 

In  general  the  controllability  gramian  V{xm)  can  be  computed 

OO 

V{xm)=  Y  Ar(x_J)m^'(x.i)m.  (1.99) 

j=— m+1 

Equation  (1.100)  can  be  further  averaged  as  in  [43]  to  include  different  initial  condition  matrices  and 
different  scalings  of  the  initial  conditions.  The  computation  of  the  empirical  controllability  gramian  is 
summarized  as  follows. 

Lemma  1.2.1  Let  the  matrix  N{x~j)m  be  the  state  responses  at  time  index  m  resulting  from  the 
collection  of  impulse  inputs  V  applied  at  time  index  — j,  where  the  £-th  column  ofAr(x-j)m  corresponds 
with  the  l-th  column  ofV.  Then  the  empirical  controllability  gramian  at  xm  is 

OO 

P(lm)  =  Y  V(x_j)mV'(x_j)m.  (1.100) 

m+1 


Note  that  for  any  j  and  m,  from  (1.95), 

>A/*(x_j) 77*4-1  i4(xj7*)y4(iT77i_i) . . .  ^4(x-.jf)^(x—^)  A{xjji)Af (x— j'jm  (1.101) 

Corollary  1.2.2  The  empirical  controllability  gramian  can  be  determined  by  the  iteratively  through 
the  Lyapunov  equation 

'Pfam+l)  ~  N (%m)m+ (^m)m-hl  +  ^(®m)^(®m)-^  (^m)  (1.102) 

Proof.  We  have  from  (1.101) 

OO 

'P(.Xm+ 1)  =  Y  ^(x-i)m+lAf'(x-j)m+i 

j=-m 

OO 

=  Af(xm±\)mAf  (^m-hl)m  A{x  771 )  E  N{X—j)m+l’hf  (x_j)m+ij4  (xm) 

j=— m-t-1 

N (^m)m-f  1  d"  A{Xrr\fP(Xrrx)A  (x^) 


(1.103) 


Of  course,  the  first  term  in  (1.103)  is  A/’(xm)m+i7\/'/(xm)m+i  =  B(xm)Bf(xTn ),  so  Corollary  1.2.2 
states  that  the  empirical  gramian  can  be  computed  by  simulating  the  previous  data  an  additional  time 
step  and  adding  the  state  responses  from  a  single-step  simulation. 

Reverse-Time  Adjoint  System 

As  a  side  note,  an  alternative  formulation  of  the  empirical  controllability  gramian  can  be  obtained  by 
simulating  the  reverse-time  adjoint  system 

7fc-i  =  A'ixkhk, 

tpk  =  B'{xk)  7fc, 

with  a  series  of  non-zero  initial  conditions,  as  in  the  empirical  observability  gramian. 

Computations  for  Empirical  Gramians 

Computation  of  the  empirical  gramians  is  similar  in  spirit  to  [43  and  [62]  with  some  key  differences. 
First,  the  formulae  in  [43]  subtract  the  temporal  means  from  the  trajectories,  which  is  not  applicable  in 
the  current  setting.  Indeed,  in  the  current  setting  we  must  track  the  trajectories  of  the  nominal  mean 
states  and  outputs.  Secondly,  the  Lemma  5  in  [43]  has  the  implicit  assumption  of  time-invariance,  as 
does  [62]. 

Perturbation  Trajectories 

Consider  the  simulation  of  the  system  (1.84)  with  a  fixed  control  input  {u*}  with  a  set  of  N  initial 
conditions  {xo}j=iAr  with  average  xq .  We  have,  for  the  nominal  trajectory  and  for  each  j  =  1, . . . ,  N 

Xk+i  =  f{xk)  +  S(^A:)>  Vk  =  h(xk)  (1.106) 

4+1  =  /(4)  Vk  =  M4)-  (1.107) 

Let  tZfc  be  a  nominal  control  input  and  let  Xk,yk  be  the  resulting  response  of  the  state  and  output, 
respectively.  We  consider  the  system  with  either  a  small  perturbation  on  the  initial  condition  t]q  or 
a  small  perturbation  on  the  control  v*  and  the  corresponding  state  and  output  perturbations  r]k  and 
<pk,  respectively.  Using  the  first  order  terms  of  a  Taylor  expansion  results  in 

Tlk+ i  =  Vk  +  ~Q~vk  Mzi OVk  +  B(uk)Vk,  (1.108) 

4>k  =  :=  C(^fc)7?fc-  (1.109) 


(1.104) 

(1.105) 
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Generalized  Gramians 

Generalized  gramians  appear  in  [9,  34]  and  are  characterized  by  the  Lyapunov  inequalities  [42] 

V(xk)  >  B(xk-\)Br(xk-\)  +  A(xk-i)V(xk-i)A'(xk- 1),  (1.110) 

Q(x*-i)  >  ^4'(arfc-i)Q(xik)i4(x/c-i)  +  C'(xk-i)C(xfc-\)*  (1.111) 

Existence  of  solutions  appears  in  [25]. 

The  empirical  generalized  gramians  can  be  computed  iteratively  through  the  time  indices:  forward 
in  time  for  the  controllability  gramians  (see  Corollary  1.2.2)  and  backward  in  time  for  the  observability 
gramians. 

The  idea  here  is  to  use  a  single  set  of  initial  conditions  and  check  that  the  matrices  of  the  resulting 
state  responses  remains  invertible,  so  it  can  be  used  for  the  computation  of  the  observability  gramian 
at  the  next  time  step.  This  way,  we  can  apply  the  Lyapunov  inequality  to  the  empirical  gramians  to 
establish  that  they  are  generalized  gramians. 

Periodic  Systems 

For  periodic  system,  the  iterative  construction  of  the  empirical  gramians  must  be  modified  due  to  the 
periodic  structure  of  the  system  matrices.  Consider  the  controllability  gramian  of  a  periodic  system 
with  period  N  so  that  V(xk)  =  V{xk+N)- 

Suppose  we  have  computed  a  nominal  empirical  gramian  V{xq)  and  we  apply  the  following  iterative 
assignments  for  k  =  0 . . .  N  —  1.  Given  V(xk)  and  6k+i  >  0  we  compute 

B(xk+ 1)  •=  B(xk)B'{xk)  +  A(xk)V(xk)A'(xk)  +  (1.112) 

>  B(xk)Bf(xk )  +  A(xk)V(xk)A'(xk)-  (1.113) 

Clearly  the  iterative  assignment  satisfies  the  Lyapunov  inequality  (1.110)  for  k  =  0 . . .  Ar  —  1  but  not 
necessarily  for  k  =  N  because  in  general  V{x^j)  V(xq)  .  We  determine  scaling  factors  m*  such  that 
V(xk)  =  m kP{xk)  satisfies  (1.110)  for  all  A;  =  0 ...  N. 

We  first  state  two  useful  facts  regarding  scaling. 

Lemma  1.2.3  If  P  >  0  and  P  >  TPT*  then  for  any  matrix  R  >  0,  there  exists  a  positive  constant 
m  such  that 

mP  >  mTPTf  +  R .  (1.114) 

Proof.  The  result  follows  by  dividing  (1.114)  by  m  and  applying  m  — ►  oo  and  P  >  TPTf.  m 

Lemma  1.2.4  If  P\  >  0  P2  >  0  and  P2  >  TPiT'  +  R  with  R  >  0,  then  ifm\  >  1  there  exists  m2  <  m\ 
such  that 

m2P2  >  miTPiT'  +  R.  (1.115) 
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Proof.  Take  m2  =  mi  to  apply  P2  >  TP\Tl  +  R  >  TPiV  +  ^R  to  obtain  the  result.  ■ 

Finally,  repeated  application  of  (1.112)  with  =  0  yields 

V{x\)  =B(xo)B'(zo)  +  A(x0)P(xo)A'(x0)  +  £1  (1.116) 

V{x  2)  =B(xi)B\xi)  +  A(xi)V(x\)Af(xi)  +  €2 

=£(xi)£'(xi)  +  A(xi)B(xo)B'(xo)A'(xi)  +  e2  +  A(xi)ei  A'(xi) 

4-  A{x\)A{xq)V{xq)A\xq)A,{x\) 


V(xn)  =£(xat-i)J3'(xjv_i)  -I - h  [A(xat-i)  . . .  A(xi)]  B(x0)B'(xo)  [A(x/v_i) . . .  A(xi)]' 

+  A(xtf-i)ejv-iA' (xjv-i)  -1 - h  (j4(xtv_i)  . . .  A(xi)]  €\  [A(xat_i) . . .  A(xi)]' 

+  [A(xyv-i)  •  •  •  ^4(^o)]  [A(x/v-i) . . .  A(xo)]'  (1-117) 

:=R  +  TV(x0)T\  (1.118) 

where 

R  :=  5(x.v-i)B'(xjv-i)  H - h  (A(x^-i)  • .  •  Mxi))  B(xo)B'(xq)  [A(x;v_i) .  * .  j4(xi)]' 

+  i4(xAT-i)e7v-i  A'(xat-i)  H - h  [A(x/v_j) . . .  A(xi)}  e\  [A(xjv-i) . . .  -A(xi)]'  (1.119) 

T  :=  [A(xjv-i)  -  *  •  ^(2:0)]  •  (1.120) 

Proposition  1.2.5  For  the  N-penodic  system ,  given  P(xo)  and  {V(xk)}k^ir..,N  specified  by  (1.110), 
there  exist  positive  constants  {mk}k=o,...,N  su°h  that  {P(xk)  :=  TrikP(xk)}k=Q%...9N-i  satisfy  (1.110)  for 
all  k  =  0, . . . ,  N  and  thus  form  the  matrix  blocks  of  a  block- diagonal  generalized  controllability  gramian 
for  the  periodic  system. 

Proof.  The  proof  is  by  construction.  If  P(xo)  >  P(xjv)  then  take  m^  =  1,  V/:  =  1, . . . ,  N  —  1.  Suppose 
V(x 0)  <  V(xn).  Since  the  nonlinear  system  has  Lyapunov  exponents  with  magnitudes  strictly  less 
than  imity,  P(xo)  >  TP(xo)Tf .  By  Lemma  1.2.3,  determine  mo  such  that 

moP(xo)  >  R  +  TmoPixcOT' .  (1.121) 

Apply  Lemma  1.2.4  to  equations  (1.116-1.117)  to  determine  {mk}k=it...tN-\  with  m^+i  <  m^  <  mo 
such  that  (1.110)  holds  for  k  =  0, . . . ,  N  —  1.  A  final  application  of  (1.121)  with  V(xn)  :=  moP(xo)  = 
'P(xo)  results  in 

P(xjv)  :=  moP(xo)  >  R  +  Tm0P(xo)T' 

=  B(xx-i)B'(xtf-i)  +  A(x7v-i)P(xjv_i)A'(xAr-i),  (1.122) 

which  is  just  (1.110)  with  k  =  N.  m 

Corollary  1.2.6  For  the  N-periodic  system ,  given  Q(xjv)  and  {Q(x^)}fc=oT...,Ar-i  specified  by  (1.111), 
there  exist  positive  constants  {qk}k=o,...,N  such  that  (Q(xa:)  :=  qkQ{xk)}k=i,...}N  satisfy  (1.111)  for  all 
k  =  0, . . . ,  N  and  thus  form  the  matrix  blocks  of  a  block-diagonal  generalized  observability  gramian  for 
the  periodic  system . 
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1.3  Stability  Analysis  of  Systems  with  Symmetry-Breaking 

1.3.1  Absolute  Stability  of  Coupled  Dissipative  Parabolic  Equations  with  Wave- 
Speed  Mistuning 

Recent  work  has  focussed  on  the  stabilizing  properties  of  symmetry-breaking  in  oscillator  systems. 
We  consider  the  problem  of  achieving  global  absolute  stability  of  an  unstable  equilibrium  solution  of 
coupled  dissipative  parabolic  equations  with  non-homogeneous  coefficients.  In  particular,  we  consider 
the  stabilization  of  a  PDE  model  describing  thermo-acoustic  instabilities  with  wave-speed  mistuning. 
Sufficient  conditions  for  absolute  stability  of  the  infinite-dimensional  system  are  established  by  the 
feasibility  of  two  finite-dimensional  linear  matrix  inequalities  (LMI).  Numerical  results  are  presented 
for  an  example  problem. 

We  establish  conditions  for  absolute  stability  of  the  nominally  unstable  zero  solution  of  a  coupled 
parabolic  equation.  In  particular  we  are  interested  in  the  stabilizing  effects  of  wave-speed  mistuning 
on  thermo-acoustic  systems  with  skew-symmetric  nonlinear  coupling.  Some  analysis  of  the  finite¬ 
dimensional  truncation  of  this  system  has  appeared  in  [33,  71,  26]  and  the  results  have  been  demon¬ 
strated  in  an  actual  engine  [17].  In  most  combustion  dynamic  applications,  identifying  the  exact  form 
of  the  heat-release  coupling  is  difficult,  and  therefore  it  is  appropriate  to  study  the  absolute  stability 
properties  of  such  systems,  where  the  skew-symmetric  feedback  is  characterized  by  a  sector  bound. 
Recent  work  has  revealed  the  importance  of  skew-symmetric  coupling  in  wave-equation  systems  [7,  8] 
and  how  mistuning  can  enhance  the  stability  of  such  systems.  This  phenomena  is  also  common  in 
suppression  of  blade  flutter  instabilities  [10,  60,  58,  70].  Additionally,  recent  work  has  focussed  on 
analysis  of  heterogeneous  distributed  systems  [24], [40]. 

Through  Lyapunov  analysis  we  establish  sufficient  conditions  for  absolute  stability  by  the  feasi¬ 
bility  of  a  set  of  finite-dimensional  linear  matrix  inequalities  (LMI)  [31,  32].  This  paper  is  organized 
as  follows.  In  section  1.3.1  we  describe  the  model  of  thermo-acoustic  instabilities  with  wave-speed 
mistuning.  This  model  serves  as  a  motivating  example  to  consider  more  general  systems  that  have 
off-diagonal  linear  coupling  which  is  described  in  section  1.3.1.  Stability  analysis  of  the  nonlinear 
infinite-dimensional  system  is  presented  in  section  1.3.1.  In  section  1.3.1,  a  finite  number  of  LMI  fea¬ 
sibility  conditions  are  presented  that  are  sufficient  for  infinite-dimensional  absolute  stability  and  the 
relation  with  finite-dimensional  model  truncation  is  presented  in  section  1.3.1.  A  numerical  example 
is  shown  in  section  1.3.1. 


Thermoacoustic  Model  with  Mistuning 


We  consider  the  following  system  that  describes  rotating  thermo-acoustic  instabilities  as  described  in 
[33,  26],  with  mistuning  parameter  a2 (9)  that  serves  to  couple  different  modes  [7,  50]. 


du 

~dt 

dp 

dt 


2,n\dp  d2 
■°  mas +  vWn 

du  d2 

-ao  +  (Wp+hu+9{u)- 


(1.123) 

(1.124) 


for  9  €  [0,  2tt)  with  periodic  boundary  conditions.  In  (1.123,1.124)  u  denotes  the  transverse  acoustic 
velocity,  p  denotes  the  acoustic  pressure,  a2 (6)  denotes  the  acoustic  wave-speed,  v  and  £  represent  vis- 
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cous  and  acoustic  damping,  respectively,  and  K  and  g{-)  denote,  respectively,  the  linear  and  nonlinear 
components  of  the  heat-release  coupling  (see  e.g.  [55,  18]).  Define  the  state  space 


x  = 


u 

P  J 


(1.125) 


It  is  convenient  to  cast  (1.123,1.124)  as  an  evolution  equation,  so  we  define  the  operators 


A  = 

so  (1.123,1.124)  become 


,,  02  „2  a 

vm  ~a  m 

—§d+K 


d 


,C=[1  0], 


— x  =  Ax  +  Bg  ( Cx ) . 


(1.126) 


(1.127) 


We  now  define  some  finite-dimensional  matrices  and  vectors  that  will  be  used  in  the  infinite¬ 
dimensional  stability  analysis.  We  will  use  {sin,  cos}  basis  functions  and  denote  the  k-th  basis  function 
as 


MO)  ••=  4= 


sin  (k&)  0  cos(A:0)  0 

0  cos  (kd)  0  sin  (kd) 


find  the  fc-th  state  vector  is 


r2n 

xk:=(Mx)-  <t>l(6)x(0) 

Jo 


d6 


(1.128) 


(1.129) 


Here  the  superscripts  [•]*,  [-]c  denote  the  sin  and  cos  components,  respectively.  As  in  [26,  50]  we 
consider  the  mistimed  wave-speed  of  the  form 


a2  (9)  =  do  +  a|  cos(20). 

For  k  =  1  the  finite-dimensional  projection  of  A  is 


(1.130) 


—  V  (Xq  —  CL2  0 


A\  :=  (<Pi,A<t>i)  = 

For  k  >  1  define  the  matrix  Ak  €  R4x4  such  that 

Afc  :=  {ipk'AM  = 


-1 

0  0 

K  0 


K 


0 

0 


~  IA  — Ctn  Go 


(1.131) 


— k2u  ka,Q  0  0 

-k  -k2Z  K  0 

0  0  —k2u  —  fcajj 

K  0  k  - k 2£ 


(1.132) 


Analysis  of  the  system  truncated  to  a  single  pair  of  modes  has  appeared  in  [33,  26] ,  where  the  mistuning 
produces  beneficial  coupling  between  the  first  pair  of  modes  [7].  The  mistuning  also  produces  coupling 
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between  mode  pairs  k  and  k  +  2.  For  k  >  1,  define 


and  for  k  >  1 


:=  -A<j>k) 


0  fca2  0  0 

0  (f  0  0 

0  0  0  —ka2 

0  0  0  0 


Afc  :=  (<j>k+2, -A<pk) 


0  ka%  0  0 

0  0  0  0 

0  0  0  —  ka2 

0  0  0  0 


Finally,  we  write  the  exogenous  input  and  output  matrices  for  the  nonlinear  feedback, 


(1.133) 


(1.134) 


and 


B  = 


'  0 

0  ■ 

1 

0 

0 

0 

.  0 

1 

1 

r2n 

J- 

Jo 

,  <?  = 


0  0  10 
10  0  0 


cos (kd) 
sin(A;0) 


g(Cx(9))  dd. 


(1.135) 


(1.136) 


Assumptions  and  Notations 

The  mistuned  thermo-acoustic  model  described  above  serves  as  a  specific  example  motivating  the 
study  of,  possibly  infinite-dimensional,  systems  with  the  following  structure. 


±i  =  A\x\  +  Ajx2  +  Bg(y)  (1.137) 

ik  =  +  Akxk  -f  A^+1Xfc+i  +  Bg(y), 

k  =  2, . . .  ,oo 
y  =  Cx 

x:=[xj  x\  ...]T.  (1.138) 

where 

Ak  =  A  +  kA W  +  k2A&  (1.139) 

A£  =  A(0l+)  +  A:A(1>+)  -I-  k2 A(2’+)  (1.140) 

A*  =  A*0-"5  -1-  A:A(1  +  fc2A<2  (1.141) 
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Referring  to  (1.132),  for  the  thermo-acoustic  model  we  have 


A  ■.= 


:= 

A<2>  := 


0  0  0  0 
0  0  K  0 
0  0  0  0 
K  0  0  0  . 

0  a2  0  0  • 

-10  0  0 
0  0  0  —Oq 

0  0  1  0  . 

-v  0  0  0 

0  0  0 

0  0  -1/  0 

0  0  0 


(1.142) 


(1.143) 


(1.144) 


Similarly,  referring  to  (1.133,1.134),  for  the  thermo-acoustic  model  we  have  A(0,+)  =  A^0,  1  = 
0,A<2'+>  =  A(2’->=0,and 


A(i-+)  _  a(j  = 


0  a2  0  0 

0  0  0  0 

0  0  0  —a2 

0  0  0  0 


(1.145) 


The  A  blocks  represent  coupling  between  modes  with  different  indices.  As  shown  in  (1.137)  the  [•]+ 
superscript  indicates  the  influence  of  mode  k  -f  1  on  mode  k  and  the  [•)“  superscript  indicates  the 
influence  of  mode  k  —  1  on  mode  k. 

The  linear  part  of  the  system  has  a  banded  infinite-dimensional  matrix  representation 


A  := 


A\  Aj  0  0 

AJ-  A2  A+  ••• 

0  AJ  A3  A  X  0 


(1.146) 


0 


0 


where  all  of  Ak, A^,A^  are  the  same  size.  In  section  1.3.1  this  restriction  is  partially  removed  by 
allowing  a  larger  size  A\  and  accordingly  zero-padding  the  matrices  Aj  and  Aj“. 

The  infinite-dimensional  matrix  representation  of  a  linear  operator  V  :  L2(Q)  x  L2(Q)  — ►  L2(Q)  x 
Z/2(fl)  consists  of  elements  Phk  G  Rnxn,  for  j,  k  =  1,2, ... ,  that  are  defined  as  follows  (see  e.g.  [1]).  If 
V  is  characterized  by  the  kernel  V{d>,  0),  so  that 

{Vx){e')=  [  v(e,.0)x(e)  de 

J  n 
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then 


:=  f  [  ${e')v{e\e)<j>k{e)de'  de. 

JaJn 


(1.147) 


For  x,  y  €  L2(Q)  we  have 

(y, Px)  =  E  E  piMxk-  (1.148) 

j  k 

Suppose  that  Pkyk  =  and  that  =  0  when  j  ^  k.  Hence  the  infinite-dimensional  matrix 
representation  of  V  is  block-diagonal,  and  we  refer  to  the  operator  V  simply  as  block-diagonal  with 
block  elements  P*.  Then  for  x,y  E  Z,2(H)  we  have 

(y,Vx)  =  J^y^PfcXfc.  (1.149) 

k 


To  show  absolute  stability,  we  will  use  a  block-diagonal  operator  P,  characterized  by  the  finite¬ 
dimensional  matrices  Pk,  in  a  Lyapunov  function.  It  is  easy  to  check  that  if  there  exists  a  constant 
a  >  0  such  that  for  all  k  =  1, 2, . . .  the  positive  definite  matrices  Pk  >  crl  then  P  is  coercive  (see  [31]). 
This  fact  will  be  applied  in  the  forthcoming  sections. 

We  first  consider  only  the  linear  system 


Wtx  =  Ax ■ 


(1.150) 


Consider  the  Lyapunov  function  V (x)  =  (x,  Px)  and  assume  that  P  has  a  block-diagonal  infinite¬ 
dimensional  matrix  representation 


P  := 


Pi  0  0 

0  P2  0 
0  0  P3 


(1.151) 


where  Pk  —  (4>k,P<Pk)  is  the  same  size  as  Ak-  Taking  the  time  derivative  of  V(x)  yields 


where 


V  ={x,VAx)  4-  (Ax,Vx) 


=^EXH PkM  +  AlPk)  Xk 

~  A:>1 


+£ 


Xk- 1 
Xk 


T 

Mk 


Xk-1 

Xk 


Mh  = 


1  °k' 

Gk  *k 


(1.152) 


(1.153) 


(1.154) 
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and 


(1.155) 

(1.156) 


Gk  :=  A+rPfc_i  4-  Pk A-_j 
:=  rnk  {PkAk  +  AkPk) 

mt  =  {  !;*='.  (1.157) 

It  is  clear  that,  under  suitable  conditions,  if  for  all  k  =  1, 2, ... , 

PkAk  +  AlPk  <  -e  (1.158) 

and 

Mk  <  0  (1.159) 

then  the  linear  system  (1.150)  is  asymptotically  stable.  We  will  discuss  the  suitable  conditions  in  more 
detail  in  the  next  section. 


Infinite-Dimensional  Absolute  Stability 

Next,  we  consider  the  system  with  nonlinear  coupling.  We  assume  that  the  nonlinear  function  g(u) 
satisfies  the  spatially  constant  sector  condition,  with  fi>  0, 


{g(u),g(u)  —  fiu)  <0,  Vu  €  R. 


(1.160) 


Note  that  the  case  of  a  spatially  varying  sector  bound  is  treated  in  [31]. 

Consider  the  Lyapunov  function  V(x)  =  (x,Vx)  ,  where  the  bounded  linear  operator  V  is  coercive 
and  block  diagonal  and  the  matrices  Pk  are  as  defined  above.  Due  to  the  coercivity  of  V,  there  exists  a 
matrix  P  >  0  such  that  Pk  >  Pyk .  Prom  (1.149),  (x,Vx)  =  Taking  the  time  derivative 

of  1/,  where  the  necessary  regularity  conditions  are  outlined  in  Chapter  6  of  [75],  results  in 


V  =  (x,VAx)  +  {Ax}  Vx)  —  2  (x,VBg(Cx)) 

<  (x,VAx)  +  (Ax,Vx) 

-  2(x,  VBg{Cx))  -f  2([fiCx  -  g{Cx)\ ,  g{Cx)) 

r  i  T  r  l 


=  E 

k>2  l 


Xk-l 

Xk 


Mk 


Xk-l 

Xk 


+  lY,xk  [PkAk  +  AlPk]  xk 


k>  1 


-  2(x,  VBg{Cx ))  +  2(\nCx  -  g{Cx)),  g{Cx)) 


=  E 


^fc-i 

xk 


Mfc 


Xk-l 

Xk 


+  ^ o  [PkAk  +  -^fc] Xk 

k>\ 

+  2x\  {vCT  -  PkB)  gk  -  2 gkgk. 


(1.161) 

(1.162) 


(1.163) 
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The  estimate  (1.162)  is  obtained  by  applying  the  sector  bound  (1.160).  For  notational  convenience 
we  define 


Vk  [PkM  +  AlPk]  xk  +  2xl  (pCj  -  PkBe )  gk 


2 

-  2 gl  9k, 


(1.164) 


so  according  to  (1.161-1.164), 


k>  2 


Sfc-1 

Xk 


Mk 


Xk-1 

Xk 


(1.165) 


It  is  now  clear  that  if  Mk  <  —s  and  14  <  0  for  all  k  then  the  system  (1.127)  will  be  asymptotically 
stable.  The  terms  Mk  appear  due  to  the  linear  coupling  from  mistiming  and  the  terms  14  appear  due 
to  nonlinear  coupling.  Each  term  contains  a  linear  stabilizing  part  that  is  used  to  make  the  each  of 
the  respective  terms  negative  definite. 

We  will  analyze  the  stability  of  the  infinite-dimensional  system  by  estimating  14  for  all  k  =  1, 2, - 

If  there  exist  constants  e  >  0  and  a  >  0  such  that  for  each  k ,  there  exists  a  matrix  P*  >  a  I  >  0  such 
that 


14  :=  Pfc-Ajt  +  A^Pk  +  e 
Lk  ~  —  PkB 

14  -f  LkLk  <  0 


then  (1.164)  becomes 

Vk<-  e^XkPkXk  -  LkLkXk  +  2xlLkgk  -  2glgk 
<  ~  e^xlPkxk  -  ~x^LkLlxk  + 

1  r 

^  ~  £^xkPkXk. 

We  can  rewrite  (1.168)  as  an  LMI  [16] 


Pk  >  (?I, 

which  we  refer  to  as  the  strictly  positive  real  LMI. 


Ft  Lk 


Ll  -/ 


<0, 


(1.166) 

(1.167) 

(1.168) 


(1.169) 


Lemma  1.3.1  If  for  k  =  1  the  LMI  (1.169)  is  feasible  and  for  k  >  2  the  LMI  (1.159)  and  (1.169) 
are  feasible,  then  system  (1.127)  is  asymptotically  stable. 

Proof.  Feasibility  of  LMI  (1.169)  for  k  >  1  implies  that  Vk  <  —e^x^xk-  For  k  >  2  feasibility  of  LMI 
(1.159)  implies  that  Mk  <  0.  Following  (1.165), 

V<-\e{x,x).  (1.170) 
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Finite  Set  of  LMI  Conditions 

For  practical  analysis,  one  would  like  to  avoid  checking  the  feasibility  of  an  infinite  number  of  LMI. 
Given  the  structure  of  the  subsystems,  it  is  possible  to  establish  sufficient  conditions  so  that  the 
feasibility  of  only  a  finite  number,  N  <  oo,  of  LMI  need  be  checked.  The  sufficient  conditions  establish 
that  if  the  LMI  are  feasible  for  k  =  N  (for  some  fixed  N )  then  they  will  be  feasible  for  all  k  >  N. 

First  fix  some  N  >  1  which  serves  as  a  truncation  index.  We  are  concerned  with  indices  k  >  N  so 

we  take  k  =  N  +  n  with  n  >  1.  Now,  for  any  n  >  1,  following  (1.139-1.141)  we  have 

AN+n  =  A  +  (N  +  +  (N  +  n)2A(2) 

=  A  +  NA[1)  +  N2A&  +  nA(1)  +  (2  nN  +  n2)A(2) 

=  Aft  +  4-  (2nN  +  n2)A^  (1.171) 

Similarly, 

Aiv+n  =  +  nA(1’+)  +  (2  nN  +  n2)A(2>+)  (1.172) 

&N+n  =  An  +  +  (2niV  4-  n2)A(2>">.  (1.173) 

Suppose  further  that  Pjv+n  =  Pn  :=  P.  The  finite  set  of  feasibility  conditions  will  be  given  for  the 
coupling  LMI  (1.154,1.159)  and  the  strictly  positive  real  LMI  (1.169). 

Coupling  LMI  (1.154,1.159) 

Define 


£(1)  ._  A(l,+)Tp  +  pA(l.-) 

(1.174) 

G(2)  ;=  A(2,+)Tp  +  pA(2,-) 

(1.175) 

tf(1)  :=  mN  +  A(1)7>) 

(1.176) 

tf<2)  :=  mN  (PA(2)  +  A(2)rp)  , 

(1.177) 
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and  note  that  mjv+n  =  mpi  =  Then, 

Gflr+n  —  A]y^nP  4-  PA^+n_j 
=  &%tP  +  PA-n_x 
+  n  (a{1-+)7>  4  PA(1’~>) 

4-  (2 nN  +  n2)  (a(2>+>tP  +  PA<2 

=  Gn  +  tiG^  +  (2  nN  4  n2)G™ 
VN+n  =  mN  (PAN+n  +  Aj,+nP) 

=  mN  (. PAn  4-  A$rP) 

+  mjvn  (PA(1)  +  v4(1)Tp) 

4  mN(2nN  4-  n2)  (PA<2)  +  vl(2>rP) 
=  VN  +  +  (2  nN  4  ri2)^2) 


and 


^+*-1  =  4  4  (2 n{N  -1)4  n2)^(2> 

=  ^at_i  4  n^(1)  -  2n^(2)  4-  (2 nN  4  n2)^(2). 


Following  (1.154), 


Ml V+n  — 

=  Mfj  4n 
4  (2niV  4  n2) 


^.v+n-i  GTN+n 

Gtf+n  ^N+n 

^(1)_2^(2)  G^t 

gw  yi*) 

^(2)  G(2)T 

G(2)  tf<2) 


It  is  easy  to  show  the  following: 
Lemma  1.3.2  Suppose  that 


3,(2)  G(2)r 


G<2) 

*(2) 

<0 

and  for  some  n  >  1, 

1 

\J/( l)  _  2^(2) 

GWr  - 

i  _1_ 

'  \p(2) 

GW  ‘ 

21V  +  n 

G(1) 

tfU) 

-T 

G<2) 

\p(2) 

Then  j  >  n  implies  that 

1 

4>(1)  _  2^(2) 

G(1)r  ' 

4 

\J<(2) 

GC2)^  ‘ 

21V  4; 

G<» 

^(i) 

G<  2) 

\J/(2) 

<  0. 


<  0. 
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(1.178) 

(1.179) 

(1.180) 

(1.181) 

(1.182) 

(1.183) 

(1.184) 


Corollary  1.3.3  Suppose  that  Mw  <  0  and  that  forn  =  1  the  LMI  (1.182,1.183)  hold.  Then  < 
0  for  all  j  >  1. 


Strictly  Positive  Real  LMI  (1.169) 
From  (1.166)  and  (1.171)  it  follows  that 


IV+n  P +  A]y+nP  4-  e 
=  IV  +  n  (P4(1>  +  AWTP^j 

+  (2 nN  +  n2)  (Pyl(2)  +  yl^p)  (1  185) 

As  in  lemma  1.3.2,  it  is  easy  to  show  the  following: 

Lemma  1.3.4  Suppose  that 

PAW  +  aV)tP<0  (1.186) 

and  for  some  n  >  1, 

— ^  (pa(1)  +  yl(1)rp)  +  (pyl(2)  +  4(2)tp)  <  0.  (1.187) 

Then  j  >  n  implies  that 

— [—  (p.4(1)  +  4(1)rp)  4-  (pA(2)  +  4(2)rp)  <  0.  (1.188) 

Corollary  1.3.5  Suppose  that  for  index  k  =  N ,  P  >  ol  satisfies  (1.169)  and  (1.186,1.187)  with 
n  —  1.  Then  P  is  a  feasible  solution  of  the  LMI  (1.169)  for  all  k  >  N  . 

Proof.  Because  Ptf+n  =  P  is  constant  over  n,  the  only  entry  in  (1.169)  that  varies  with  n  is  the  (1, 1) 
entry.  Feasibility  of  (1.186,1.187)  ensures  that  IV+j  <  T^v-fn  when  j  >  n.  m 


Model  Truncation 


Restricting  the  Lyapunov  fimction  to  be  block  diagonal  can  be  overly  conservative.  This  restriction 
can  be  partially  reduced  by  allowing  allowing  larger  block  sizes.  In  this  section  we  consider  a  larger 
finite-dimensional  model  and  analyze  the  feasibility  of  the  LMI.  This  allows  at  least  the  first  block 
to  be  arbitrarily  large,  but  finite.  Extending  the  analysis  to  apply  to  larger  subsequent  blocks  is 
si  raightforward  and  not  explicitly  discussed  here. 

For  some  N  >  1  define  truncated  finite-dimensional  model  Ab  €  R(iV-1)nx(^~1)nj 


Ai 

A+ 

0 

0 

A*i 

a3+ 

0 

^2 

a3  a+ 

0 

*N-l 

0 

O 

> 

1 

to 

An- i 

(1.189) 
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and 


Bb  —  diag{B, . . . ,  B },  Cb  —  diag{C, . . . ,  C } 


(1.190) 


xD:=[x{  ...  ]T  (1.191) 

5s  :=  [  5iT  •  •  •  9ji- 1  ]  •  (1.192) 

In  each  of  (1.189,1-190)  the  0  blocks  are  assumed  to  be  appropriately  sized  matrices.  Consider  a 
block-diagonal  Lyapunov  function 

V'(x)  =  (x,Vx)  =  xtbPbxb  +  xkpkxk  (1.193) 

k>K 


where  PD  E  and  Pk  €  Rnxn  for  k  >  N.  Following  (1.154-1.156)  define 


where 


and 


r*B  gtd  1 

B  [Gb  J 

(1.194) 

Gb  :=  A  +TPB  +  Pn&b 

(1.195) 

*b-=\{PbAb  +  AIpb) 

(1.196) 

:=  [  0  ...  0  A^r  ]  €  Rnx(N_1)n 

(1.197) 

Ab  :=  [  0  ...  0  A^_!  ]  €R»*(*-i)n 

(1.198) 

and  \  {Pn-An  +  Aj,PN)  by  (1.156). 

Taking  the  time  derivative  of  V{x)  and  applying  the  sector  bound  (1.160)  results  in 


V  <  xb~  [PbAb  +  A'bPb\  xb 


+  2 xf  (gCB  -  PdBb)  gB  -  2 gB9D 

T 


+ 


Xjg 

xn 


+  E 

k>N+l  L 


Mb 

Xk 


Xb 
XN 

iT 

Mh 


Xk—1 

Xk 


+  £**5  [pkAk  +  AlPk]  xk 


k>N 


+  2xl  {vCT  -  PkB )  gk  -  2g{gk. 


(1.199) 

(1.200) 

(1.201) 

(1.202) 

(1.203) 

(1.204) 


Define 


Tb  :=  PqAb  +  AqPb  +  e 
Lb  ■=  nCg  —  PbBb 

and  consider  the  LMI  analogous  to  (1.169) 


Pb  >  <rl, 


Tb  Lb 
Ll  -I 


<0. 


(1.205) 

(1.206) 


(1.207) 


Theorem  1.3.6  Given  some  positive  truncation  index  N  >  1.  Suppose  the  following  conditions  hold: 

1.  Pb  satisfies  LMI  (1.207) 

2.  Pn  satisfies  LMI  Mb  <  0 

3.  Pn  satisfies  LMI  (1.182,1.183)  with  n  =  1 
4-  Pn  satisfies  LMI  (1.186,1.187)  with  n  =  1 

Then  the  system  (1.137)  is  absolutely  stable  with  respect  to  feedback  nonlinearities  satisfying  (1.160). 

Proof.  Take  the  Lyapunov  function 

V(x )  =  (x,Vx)  =  x^Pbxb  +  Bn xk-  (1.208) 

k>N 


The  stated  conditions  ensure  that  the  terms  in  (1.199-1.204)  arc  negative: 

1.  implies  that  (1.199,1-200)  <  —  cx^ib- 

2.  implies  that  (1.201)  <0. 

3.  implies  that  (1.202)  <  0. 

4.  implies  that  (1.203,1.204)  <  —<J2k>Nxkxk- 

Together  these  imply  that  V  <  —  e||x||.  Furthermore,  since  the  infinite-dimensional  matrix  represen¬ 
tation  of  V  is  block  diagonal  consisting  of  the  two  finite-dimensional  blocks  Pb  and  Pn,  the  operator 
V  is  coercive  and  bounded.  ■ 
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Numerical  Results 


In  this  section  we  present  numerical  results  concerning  the  solutions  of  the  LMI  given  in  the  previous 
sections.  We  consider  the  PDE  system  (1.123,1-124)  with  wave-speed  mistuning  and  a  feedback  non¬ 
linearity  satisfying  a  sector  bound.  As  discussed  in  {7], [26]  it  is  known  that  wave-speed  mistuning  acts 
to  stabilize  thermo-acoustic  systems  with  skew-symmetric  heat-release  coupling.  However,  in  these 
cases  the  results  were  established  only  for  finite-dimensional  models.  The  model  parameters  were 
chosen  such  that  the  system  without  mistuning  is  unstable  and  develops  a  limit-cycle  oscillation.  A 
wave- speed  mistuning  parameter  along  with  truncation  index  N  was  chosen  such  that  the  conditions 
of  theorem  1.3.6  were  satisfied.  The  response  of  the  acoustic  pressure  p{0,t)  is  shown  in  figure  1.18. 
Initially  the  system  has  no  mistuning  and  develops  a  limit-cycle  due  to  the  skew-symmetric  feedback. 
At  t  =  25  the  wave-speed  mistuning  is  applied  and  the  system  stabilizes  to  the  uniform  steady  solution. 


Figure  1.18:  wave-speed  mistuning  initially  turned  off.  Turned  on  at  time=  25. 


1.3.2  Coupling  of  Stable  Subsystems 

In  many  applications,  feedback  interconnections  between  stable  dynamical  subsystems  can  lead  to 
instability  or  high  sensitivity.  Examples  of  such  instabilities  include  beam  vibration  [6,  12],  flutter 
instabilities  [?],  and  large  platoons  of  vehicles  [?,  ?,  ?,  39,  8].  Often,  the  dynamic  models  of  these 
systems  are  characterized  by  a  large  number  of  similar  subsystems  through  neighboring  coupling,  and 
the  coupling  of  the  subsystems  may  have  aperiodic  or  random  patterns  [70,  60].  As  a  result,  stability 
analysis  of  such  dynamical  systems  becomes  a  challenge. 

Under  this  contract,  we  considered  coupling  of  arbitrarily  many  stable  subsystems  connected  in  a 
chain  form  (see  Figure  1.19),  and  we  studied  how  stability  and  robustness  evolve  for  varying  coupling 
gains  between  neighboring  subsystems.  In  particular,  we  focused  on  linear  stable  subsystems  that  have 
a  property  of  “ negative  imaginary  frequency  response f  [44],  which  basically  means  that  the  subsystem’s 
phase  shift  angle  is  within  the  range  of  [— 7r,  0]  and  closely  relates  to  the  property  of  counterclockwise 
input-output  dynamics  studied  in  [2].  We  call  them  stable  negative  imaginary  (SNI)  systems.  It  has 
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been  shown  in  [44]  that  stability  for  a  positive  feedback  interconnection  of  two  SNI  subsystem s  can  be 
checked  by  their  DC  gain  (i.e.  the  loop  gain  at  zero  frequency);  see  Figure  1.20. 


Figure  1.19:  String  of  n  coupled  SNI 

The  string  of  n  SNI  subsystems  is  simply  the  coupling  of  n  SNI  subsystems  via  neighboring  pos¬ 
itive  feedback  interconnection.  Following  [44,'  2],  we  characterized  the  stability  of  the  string  through 
convergence  of  a  continued  fraction  that  denotes  DC  gain  of  sequentially  coupled  subsystems.  The 
derived  results  avoid  explicit  computation  of  eigenvalues  or  construction  of  Lyapunov  functions,  and 
they  are  readily  applicable  to  analyzing  stability  and  robustness  of  dynamical  systems  with  different 
coupling  patterns. 

Similar  analysis  was  also  extended  to  a  ring  form  (See  Figure  1.21). 
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Figure  1.20:  Positive  feedback  interconnection  of  two  subsystems  G{s)  and  $(s). 


Figure  1.21:  Ring  of  n  coupled  SNI 


51 


The  String  Form  of  Coupling 


Denote  by  $i(s)  the  transfer  function  of  the  first  (i  + 1)  coupled  subsystems  (see  the  transfer  functions 
contoured  by  dashed  squares  in  Figure  1.19).  Inductively,  we  obtain  that, 


<M*)  = 

*i(s)  = 


G2(s) 

1  -  CiriGi(s)G2(s)’ 
_ C»+i(s) _ 

1  -Cir,$i_i(s)Gi+i(s) 


Vi  =  2,3,  •  • 


,n  —  2. 


Thus  the  string  of  n  coupled  stable  subsystems  is  exactly  the  positive  feedback  interconnection  of 
$n-2(s)  sind  Cn_irn_iGn(s).  Also,  denote  the  DC  loop  gains  of  coupled  subsystems  as  follows: 


Ai  =  ciriGi(0)G2(0) 


and,  for  each  i  =  2, 3,  •  •  •  ,  n  —  1, 


Ai  :=Ciri$i^x(0)Gi+1(0) 

_  CiriGi(0)Gi+i(0) 

“l  -  Ci^n^i^GiiO) 

_CiriGi(0)Gi+i(0) 

1  -  * 

Clearly,  each  Ai  is  a  continued  fraction.  Thus  we  have 


CirjGt(0)Gi+1(0) 

^  ^  Ci— l^t— iGi— i(0)Gj(0) 

1  4-  ~  ci-2rj~2C:i^2(0)Gri^i(0) 


Vi  =  2,3,  •  •  •  ,n  —  1. 


(1.209) 


The  following  result  extends  [44,  Theorem  5]  to  the  setting  of  a  series  of  coupled  SN1  by  relating 
stability  of  the  string  to  convergence  of  the  continued  fractions. 

Theorem  1.3.7  For  the  string  of  n  coupled  SNI,  if  G\  is  a  SNI  subsystem  in  a  strict  sense , 


•  Gi(oo)  =  0  for  all  i  =  1, 2,  •  •  •  ,  n, 


•  At-  in  (1.209)  is  less  than  unity  for  all  i  =  1, 2,  •  •  •  ,  n  —  1, 
then  the  string  is  internally  stable. 

The  next  result  presents  a  sufficient  condition  in  term  of  the  subsystems'  DC  gains  and  their  local 
coupling  gains  for  guaranteeing  that  each  continued  fraction  A,*  <  1  in  Theorem  1.3.7. 

Theorem  1.3.8  The  string  of  n  coupled  SNI  is  internally  stable  if  G\  is  a  SNI  subsystem  in  a  strict 
sense ,  Gt(oo)  =  0  for  all  i  =  1,2,---  ,n7  and 

(  cinGi(0)G2(0)  <  i, 

<  anGiWGi+m  <  i  Vi  =  2,  •  •  •  ,  n  —  2,  (1.210) 

i  Cn— 1^71—1^71— 1  (0)G"n(0)  <  ^  (if  Tl  >  2). 
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The  Ring  Form  of  Coupling 

It  turns  out  that  Theorem  1.3.8  can  be  extended  to  a  ring  form  (See  Figure  1.21). 


Theorem  1.3.9  The  ring  of  n  coupled  SNI  subsystems  is  internally  stable  if  G\  is  a  SNI  subsystem 
in  a  strict  sense,  n^=i °i  =  17?=! r»>  an d  G»(0)  >  0  and  Gi( oo)  =  0  for  alii  —  1, 2, ■  •  •  , n,  and 


|  CiriGj(0)G<+i(0)  <  \  Vi  =  1,  •  •  •  ,n  —  1, 
1  c„rnGi(0)Gn(0)  < 


Robust  Stability  for  the  String 

The  quantity  An_i  as  defined  in  (1.209)  can  be  regarded  as  “robustness”  of  stability  for  the  string  of 
n  coupled  stable  subsystems.  It  is  interesting  to  see  how  varying  coupling  gains  make  An_i  approach 
its  limits  (typically  less  than  unity)  as  the  dimension  n  increases  to  infinity.  For  A„_i  <  1,  the  string 
remains  stable  with  respect  to  small  variations  of  neighbor-coupling  gains  and  certain  unmodeled 
dynamics  that  also  belongs  to  the  class  of  SNI. 

In  Theorem  1.3.8,  consider  the  case  of  c,rjGi(0)Gi+i(0)  =  for  all  i  =  1,2.  •  •  •  ,n  —  1,  where 
[i  €  [0, 1)  is  called  a  mistuning  parameter.  Solving  the  equation  A  =  •  yrj  gives  A  =  We 

verify  that  if  0  <  A<  <  then  Aj+i  >  A,-  and  lim  Ai  =  Also,  we  verify  that  if  Aj  =  —  e, 

1— ‘OO  2 

where  0  <  e  <  ^p,  then 

x  1  —  n  (1  —  ^  1  —  fx  oe 

i+1  ~  ~2  1  +  /x  +  2e  ~  2  1  +  2e’ 

where  a  =  and  the  equality  holds  for  /j.  =  0.  Consequently,  we  establish  by  induction  that,  for 
each  i  =  1,2,---  ,n  —  1, 

•  if  n  =  0  then 

\.  —  1  g 
‘  2  1  -I-  2e(i  -  1) 

and  hence  Aj  approaches  its  limit  5  in  the  order  of  <!?(i)  as  i  increases,  and 

•  if  0  <  p  <  1  then 


^  <  1  ~  M _ _ <  1  —  M 

'■  2  2  i+A 

and  hence  Ai  approaches  its  limit  in  the  order  of  at  most  O{ox)  as  i  increases. 

For  the  general  case  of  CjriGi(0)(?i+i(0)  <  (i  =  1,2, •  •  •  ,n  —  1),  the  above  convergence 

statements  still  hold,  which  is  used  for  estimating  the  values  of  the  continued  fractions  and  their 
approximants,  allows  the  parameters  |Aq|  <  1. 


Application:  Decentralized  Control  of  Vehicle  Platoons 

As  an  application,  we  considered  a  platoon  of  n  vehicles  modeled  as  follows 

yi^=Ui  Vi  =  1,2,  •  *  •  ,n,  (1.212) 

where  yi  and  m  are  the  position  and  control  input  of  the  ith  vehicle,  respectively.  The  control  objective 
is  to  maintain  a  desired  distance  L  >  0  between  any  two  neighboring  vehicles  and  ensure  a  desired 
velocity  V  >  0  for  each  vehicle,  provided  that  each  vehicle  has  the  spacing  and  velocity  information 
relative  to  its  front  and  back  ones.  Denote  6i(t)  =  j/i(£)  —  y,_i(t)  +  L. 

Two  control  strategies  are  provided  as  follows.  The  first  one  is  designed  with  control  inputs: 

Mt)  =  -  ft(yi(t)  -  V)  -  Fi{yi{t)  -  Vt) 

Ui{t)  =  “  & (yi(*)  -  V)  ~  Fi6i(t)  4-  £i<5i+i(t)  Vz  =  2,3,  •  •  •  ,n  -  1, 

^n(^)  =  £n(2/n(0  V)  —  F n^n(^)? 

where  each  &  is  a  damping  parameter  and  each  F*  (respectively,  B*)  is  a  control  gain  due  to  the 
difference  between  the  positions  of  the  zth  vehicle  and  its  front  (respectively,  back)  one.  This  control 
strategy  requires  all  vehicles  to  have  a  priori  knowledge  of  the  desired  velocity  V. 

Using  the  previous  stability  results,  we  established  that  the  system  (1.212)  achieves  asymptotic 
tracking  with  positive  real  constants  &,  F*,  and  B{  (z  =  1, 2,  •  •  •  ,  n)  is  internally  stable. 

As  the  second  control  strategy,  we  apply  the  following  control  inputs: 

ui(t)  =  -  Fi(yi{t)  - Vt  +  L )  -  kFi(yi(t)  -  V)  +  BiS2(t)  +  kBx &2(t), 
ui(t)  =  -  Fi8i(t)  -  KFi5i(t)  4-  Bi8i+i(t)  4*  KBi6i+i(t)  Vz  =  2, 3,  •  •  •  ,  n  -  1, 

=  “  Pn^n(^)  —  &F n<jn(^) 


where  k  >  k  >  0  are  constants  and  each  F*  (respectively,  Bi)  is  a  control  gain  due  to  the  spacing  and 
velocity  differences  between  the  zth  vehicle  and  its  front  (respectively,  back)  one.  This  control  strategy 
requires  only  the  leading  vehicle  has  the  information  of  the  desired  velocity  V  for  the  platoon. 

Again,  using  the  previous  stability  results,  we  established  that  the  system  (1.212)  achieves  asymp¬ 
totic  tracking  with  positive  real  constants  k>  /t,  F*  and  B*  (i  =  1, 2,  •  •  •  ,  n). 

These  stability  results  for  large  vehicle  platoon  can  serve  as  a  foundation  for  designing  various 
coupling  patterns  and/or  gains  so  as  to  enhance  robust  stability  or  improve  system  performance. 
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Chapter  2 

Transitions  of  past  AFOSR-sponsored 
research  at  UTRC 

2.1  Nonlinear  Dynamic  Texture  Modeling 

The  concepts  developed  for  nonlinear  dynamic  texture  modeling  were  applied  in  an  internally  funded 
project  to  calibrate  dynamic  models  from  synchronized  high-speed  video  and  acoustic  press.  POC 
Tory  Brogan,  Pratt  &  Whitney,  860  557  0547. 


Chapter  3 

Personnel  Supported 


UTRC  personnel:  Gregory  Hagen,  Chaohong  Cai,  George  Mathew,  Andrzej  Banaszuk,  Alberto  Sper- 
anzon  . 
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Chapter  4 

Publications 


4.1  Journal  papers 

Oi]  c.  Cai  and  G.  Hagen,  Stability  analysis  for  a  string  of  coupled  stable  subsystems  with  negative 
imaginary  frequency  response.  To  appear  in:  IEEE  Transactions  on  Automatic  Control,  2010. 

[j2)  G.  Hagen,  Stochastic  Averaging  for  Identification  of  Feedback  Nonlinearities  in  Thermo-Acoustic 
Systems.  Provisionally  accepted:  ASME  Journal  on  Dynamic  Systems,  Measurement,  and  Control, 
2010. 


4.2  Journal  papers  in  preparation 

[pi]  U.  Vaidya  and  G.  Hagen,  Model  Reduction  of  Nonlinear  Systems,  Tangent  Space  Approach.  In 
preparation 

[p2]  C.  Cai  and  G.  Hagen,  Stability  analysis  for  a  Ring  of  Coupled  Negative  Imaginary  Systems. 
In  preparation 

[p3]  G.  Hagen  and  A.  Speranzon,  Locally  Affine  Switched  Dynamic  Texture  Models.  In  prepara¬ 
tion 


[p4]  G.  Hagen,  On  Dynamic  Mode  Decomposition  and  Dynamic  Texture  Models.  In  preparation 

[p5]  G.  Hagen  and  A.  Speranzon,  Tangent  Space  Approximation  Approach  to  Compressed  Sensing 
and  Dictionary  Optimization.  In  preparation 
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4.3  Conference  papers 

[cl]  J.  Cohen,  G.  Hagen,  A.  Banaszuk,  S.  Becz,  P.  Mehta,  ’’Attenuation  Of  Gas  Turbine  Combustor 
Pressure  Oscillations  Using  Symmetry  Breaking”,  AIAA  Aerospace  Sciences  Meeting,  2011,  Orlando, 
FL. 

[c2]  G.  Hagen,  Absolute  Stability  of  a  Dissipative  Wave  Equation  with  Waves-Speed  Mistuning.  Pro¬ 
ceedings  of  the  American  Control  Conference,  Baltimore,  MD.  2010. 

[c3]  U.  Vaidya,  G.  Hagen,  Model  reduction  of  nonlinear  systems:  Tangent  space  approach.  Pro¬ 
ceedings  of  the  American  Control  Conference,  Baltimore,  MD.  2010. 

[c4]  C.  Cai,  G.  Hagen,  Coupling  of  Stable  Systems  with  Counterclockwise  Input-output  Dynamics. 
Proceedings  of  the  American  Control  Conference,  Baltimore,  MD.  2010. 

[c5]  C.  Cai,  G.  Hagen,  Stability  Results  for  String  of  Stable  Subsystems  with  Applications  to  De¬ 
centralized  Control  of  Large  Vehicle  Platoons,  Proceedings  of  the  4Sth  IEEE  Conference  on  Decision 
and  Control/  28th  Chinese  Control  Conference,  Shanghai,  China,  2009 

[c6]  M.  Arienti,  M.  C.  Soteriou,  G.  Hagen,  M.  L.  Corn,  Analysis  of  Liquid  Jet  Atomization  Dy¬ 
namics  Using  Proper  Orthogonal  Decomposition,  47th  AIAA  Aerospace  Sciences  Meeting  Including 
The  New  Horizons  Forum  and  Aerospace  Exposition,  5  -  8, Orlando,  Florida,  January  2009. 

[c7]  M.  Arienti  ,  M.  Corn,  G.  S.  Hagen,  R.  K.  Madabhushi  and  M.  C.  Soteriou,  Proper  Orthogo¬ 
nal  Decomposition  Applied  to  Liquid  Jet  Dynamics,  ILASS  Americas,  21st  Annual  Conference  on 
Liquid  Atomization  and  Spray  Systems,  Orlando,  Florida,  May  2008. 

[c8]  S.  Meyn,  G.  Hagen,  G.  Mathew,  A.  Banaszuk,  On  Complex  Spectra  and  Metastability  of  Markov 
Models,  Proceedings  of  the  47th  IEEE  Conference  on  Decision  and  Control,  Cancun,  Mexico,  2008 

[c9]  G.  Hagen,  U.  Vaidya,  An  Approach  for  Nonlinear  Model  Extraction  from  Time-Series  Data. 
Proceedings  of  the  American  Control  Conference,  Seattle,  WA.  2008. 
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